Convergence properties of a geometric mesh smoothing algorithm
CCONVERGENCE PROPERTIES OF A GEOMETRIC MESHSMOOTHING ALGORITHM
DIMITRIS VARTZIOTIS AND DORIS BOHNET
Abstract.
We describe a simple geometric transformation of triangles whichleads to an efficient and effective algorithm to smooth triangle and tetrahedralmeshes. Our focus lies on the convergence properties of this algorithm: weprove the effectivity for some planar triangle meshes and further introducedynamical methods to study the dynamics of the algorithm which may beused for any kind of algorithm based on a geometric transformation. Introduction
Preliminary remarks.
The finite element method is the standard instru-ment to simulate the behavior of solid bodies or fluids in engineering and physics.The first preparatory step of this method is the discretization of the underlyingdomain into finitely many elements which could be easily described by parameters,i.e. surfaces are mostly approximated by triangles, quatrilaterals or parallelograms.Because of the design process in modern engineering, an initial mesh for the do-main is often given, and the next important step is the preprocessing of this meshto obtain a good base for the application of the finite element method. As therequirements for simulation results are more and more strict and real time simu-lations and simulations on evolving objects present new challenges, a fast, reliableand preferably automatic mesh preprocessing is an important link in the simulationprocess.Not surprisingly, there is a wide variety of methods at hand to improve the meshquality of a given mesh, see e.g. the surveys [7] and [2]. One can identify two mainapproaches:
Geometry-based:
A geometric smoothing method changes directly the ge-ometry of the mesh, that is, it relocates the nodes. A popular example is theLaplacian smoothing which maps every node to the arithmetic mean of itsneighboring nodes (see e.g. [8] or [3]). There are also methods that changethe topology of the mesh by deleting small elements or subdividing largeones. All these methods have in common that they are usually quickly im-plemented and very fast; additionally, they can be often straightforwardlycombined with techniques of parallel computing. The main disadvantage istheir heuristic character so that the convergence of an algorithm is mostlyonly empirically, but not theoretically assured. Consequently, they aresometimes combined with optimizational approaches as in the early work[5].
Optimization-based:
The principal idea behind any optimizational approachis to define a function on the set of meshes which represents the quality ofa mesh, and to find the maximum of this function by usual numerical op-timization, e.g. gradient methods (see [6] and following articles by theseauthors). The main advantage of these methods is clearly that they lead toa mesh of higher quality, but in the case of a non-convex quality function,it can usually not be assured that the transformed mesh corresponds to the a r X i v : . [ m a t h . NA ] N ov DIMITRIS VARTZIOTIS AND DORIS BOHNET global and not only local optimum. Also, the computation is usually costlywith regard to runtime and storage. On the other hand, new advances forfast and robust solutions of optimization problems could be used as e.g.evolutionary optimization algorithms (see e.g. [16], [9]).In this article we present a geometric approach for mesh smoothing which consistsof a simple geometric transformation of every element of a mesh which do not affectthe topology. Here we only consider triangle and tetrahedral meshes. This ansatz issimilar to the GETMe algorithm introduced for triangle and tetrahedral meshes in[12, 15] and proved to be element-wise effective in [14]. But while the serial of thesearticles mainly focus on the numerical results and the improvement of runtime andperformance by adjusting the algorithm, we study the mathematics underlying ourpresented geometric method and prove that any not too distorted planar mesh oftriangles converges to the best possible mesh for the given mesh topology, that is,the difference between the normalized distances from the vertices to the centroidfor all triangles is the smallest possible. This point - the mathematical discussionof convergence properties and the application of dynamical methods - is surelythe main achievement of the present article. For completeness and as motivationfor a future application, we shortly discuss the performance of our method as asmoothing algorithm, but we do not explore the practical aspects of our algorithmin detail.1.1.1.
Organization of the article.
In Section 2 we describe the discrete geometrictransformation of a triangle element on which the smoothing algorithm is basedand discuss its mathematical properties. In the forthcoming Section 3 we derivethe smoothing algorithm from this transformation for a triangle mesh and proveits convergence for some particular triangle meshes, i.e. if the transformation isiteratively applied, the mesh converges to the best possible mesh for a given meshtopology. At the end – in Section 4 – we briefly discuss the numerical results.2.
The geometric triangle transformation
Before we start with the rigorous mathematical description we motivate thegeometric transformation, the subject of this article, and its regularizing mechanismby the following observations which have their offspring in [11]:2.1.
Introductory observations.
Imitating the rotational symmetry group action of the triangle.
The symme-try group of a regular triangle ∆ = ( z , z , z ), z i ∈ C , is the dihedral group D which is generated by a reflection and a rotation by π around the circumcenter c of the triangle. Consider the rotation: if the circumcenter lies in the origin, therotational element then acts on the triangle by mapping the vector z i − onto thevector z i for i ∈ Z .If the triangle is not equilateral, we can take the centroid, that is, the arithmeticmean of the three nodes, instead of the circumcenter and imitate the rotation bystill mapping the vector z i − onto | z i − || z i | z i such that the resulting vector has lengthequal to z i − but points into the direction of z i , i.e. we rotate the vector z i − around c .Sure, this rotation around the centroid is neither 3-periodic nor isometric, butthis action, if iterated, converges to the classical rotation by π because the centroidconverges to the circumcenter and the distances from the centroid to the verticesbecome equal. This will be shown in Subsection 2.2.2 below. ONVERGENCE PROPERTIES OF A GEOMETRIC MESH SMOOTHING ALGORITHM 3 z (1)0 z (1)1 z (1)2 z (0)1 z (0)2 z (0)0 c (0) z (1)0 z (1)1 z (1)2 c (1) z (2)0 z (2)1 z (2)2 z (2)0 z (2)1 z (2)2 c (2) Figure 1.
The first three iterations of the geometric elementtransformation represented as rotations: observe how the circleradii approach and the centroid moves to the circumcenter.2.1.2.
First intuitive explanation of the mechanism.
Consider the centroid c . Westart at the origin c (0) = 0. After the first iteration we get c (1) = (cid:16) | z || z | z + | z || z | z + | z || z | z (cid:17) . Just change it a bit setting r i := | z i − || z i | to obtain ˜ c (1) = r + r + r ( r z + r z + r z ).Observe that by the inequality of the geometric and arithmetic mean we always have3 = 3( √ r r r ) ≤ r + r + r . The new ˜ c is a weighted arithmetic mean of thevertices where the weight r i := | z i − || z i | is greater the smaller the distance | z i | com-pared to | z i − | . Accordingly, c is moved more than the average into the directionof z i if z i was much closer to c then z i − . The new distance from c (1) to | z i − || z i | z i isthen less than | z i − | . On the other hand, if z i was far from c compared to z i − , c ismoved less than the average into the direction of z i , and the new distance from c (1) to | z i − || z i | z i is greater than | z i − | . Consequently, due to the controlling weights themaximum distance from a vertex to the centroid lessens while the minimal distanceaugments such that the distances become equal in the long run. In other words,the centroid moves to the circumcenter of the triangle. It is obvious that c stopsmoving if and only if the distance to each vertex is equal or equivalently, if andonly if c is the circumcenter.In Section 2.2.2 we give a rigorous proof that this transformation iteratively appliedto any non-degenerate triangle makes it equilateral.2.2. The geometric element transformation.
Description of the geometric element transformation.
Let ∆ = ( x , x , x )with x i ∈ R or x i ∈ R for i = 0 , , c = ( x + x + x ) thecentroid of the triangle. The transformation works then as following:∆ new = ( x ,new , x ,new , x ,new ) ,x i,new = (cid:107) x i − − c (cid:107) (cid:107) x i − c (cid:107) ( x i − c ) + c for i ∈ Z . To keep the centroid fixed we move the centroid c new of the transformed triangle∆ new back into the old centroid c : x i,new = x i,new − c new + c, i = 0 , , . Combiningthese two steps into one we get for i ∈ Z with r i := (cid:107) x i − − c (cid:107) (cid:107) x i − c (cid:107) − :(1) x i,new = 23 r i ( x i − c ) − r i +1 ( x i +1 − c ) − r i − ( x i − − c ) + c. DIMITRIS VARTZIOTIS AND DORIS BOHNET
This is the whole, very simple geometric transformation which can be directlyimplemented into any mathematical software, e.g. Matlab.2.2.2.
Formal proof of element-wise convergence.
For the proof that any non-degeneratetriangle converges under Transformation (1) to a equilateral triangle we can supp-pose that the triangle lies in the euclidean plane R which we identify – for sim-plifying notations – with C : let ∆ = ( z (0)0 , z (0)1 , z (0)2 ) with z (0) i ∈ C be an arbitrarytriangle with z (0) i (cid:54) = z (0) j for i (cid:54) = j . Denote by c = ( z + z + z ) the centroid ofthe triangle. Without loss of generality, we assume that the centroid c lies in theorigin, that is, c = 0. For n ∈ N denote by r ( n ) i the ratio | z ( n ) i − || z ( n ) i | with i ∈ Z . ThenTransformation (1) becomes the following transformation, recursively defined for n ≥ i ∈ Z :(2) z ( n ) i = 23 r ( n − i z ( n − i − r ( n − i +1 z ( n − i +1 − r ( n − i +2 z ( n − i +2 We prove that the ratio r ( n ) i for i = 0 , , Theorem 1.
With the notations above, we have lim n →∞ r ( n ) i = 1 for i = 0 , , ,i.e. the ratio of the distances from the vertices to the centroid converges to . Before we start with the proof of Theorem 1 we show in the next two preliminarylemmata that the maximal distance max i | z ( n ) i | from a vertex z ( n ) i to the centroid c is a strictly decreasing sequence, and symmetrically, that the minimal distancefrom a vertex to the centroid is a strictly increasing sequence. Lemma 1.1.
For n ≥ we have max i =0 | z ( n +1) i | < max i =0 | z ( n ) i | .Proof of Lemma 1.1. We prove this Lemma by a simple estimation. Without lossof generality we assume that max i =0 | z ( n +1) i | = | z ( n +1)0 | . We have | z ( n +1)0 | = | r ( n )0 z ( n )0 − r ( n )1 z ( n )1 − r ( n )2 z ( n )2 | utilizing z ( n )0 + z ( n )1 + z ( n )2 = 0= | − r ( n )0 ( z ( n )1 + z ( n )2 ) −
13 ( r ( n )1 z ( n )1 + r ( n )2 z ( n )2 ) | < max i r ( n ) i | z ( n )0 | If r ( n )0 = max i r ( n ) i the proof is easily finished by | z ( n +1)0 | < r | z ( n )0 | = | z ( n )2 | ≤ max i | z ( n ) i | . Otherwise, assume that r ( n )1 is maximal (the case that r ( n )2 is maximal works anal-ogously). Then we substitute in the equation above z ( n )1 = − z ( n )0 − z ( n )2 and weget: | z ( n +1)0 | = | ( 23 r ( n )0 + 13 r ( n )1 ) z ( n )0 + ( 13 r ( n )1 − r ( n )2 ) z ( n )2 | r ( n )1 = max r ( n ) i .< r ( n )1 | z ( n )1 | = | z ( n )0 | ≤ max | z ( n ) i | . (cid:3) Now we prove in an analogous way that the sequence of minima is monotonicallyincreasing:
Lemma 1.2.
For n ≥ we have min i =0 | z ( n ) i | < min i =0 | z ( n +1) i | . ONVERGENCE PROPERTIES OF A GEOMETRIC MESH SMOOTHING ALGORITHM 5
Proof of Lemma 1.2.
Assume without loss of generality that | z ( n +1)0 | = min | z ( n +1) i | .We have the following estimate: | z ( n +1)0 | = | r ( n )0 z ( n )0 − (cid:16) r ( n )1 z ( n )1 + r ( n )2 z ( n )2 (cid:17) | substituting z ( n )0 = − z ( n )1 − z ( n )2 = | − r ( n )0 (cid:16) z ( n )1 + z ( n )2 (cid:17) − (cid:16) r ( n )1 z ( n )1 + r ( n )2 z ( n )2 (cid:17) | > (cid:18) r ( n )0 + 13 min i =1 , r ( n ) i (cid:19) | z ( n )0 | (3)Consequently, we have to show that (3) is greater than min i | z ( n ) i | : If r ( n )0 is minimal,we easily get (cid:18) r ( n )0 + 13 min i =1 , r ( n ) i (cid:19) ≥ r ( n )0 ⇒ | z ( n +1)0 | > r | z ( n )0 | = | z ( n )2 | ≥ min | z ( n ) i | . Otherwise, if r ( n )0 is maximal, we certainly have – utilizing the inequality of arith-metic and geometric mean (cid:18) r ( n )0 + 13 min i =1 , r ( n ) i (cid:19) ≥ (cid:18) r ( n )0 + 13 r ( n )1 + 13 r ( n )2 (cid:19) ≥ (cid:113) r ( n )0 r ( n )1 r ( n )2 = 1 , finishing the proof for this case due to | z ( n +1)0 | > | z ( n )0 | . In the last case, if r ( n )0 isneither maximal nor minimal, either r ( n )1 or r ( n )2 is maximal. Assume without lossof generality that r ( n )1 is maximal: | z ( n +1)0 | > (cid:18) r ( n )0 + 13 r ( n )2 (cid:19) | z ( n )0 | ≥ (cid:113) r ( n )0 r ( n )0 r ( n )2 | z ( n )0 | = (cid:113) r ( n )1 r ( n )1 r ( n )0 | z ( n )1 | utilizing | z ( n )0 | = r ( n )1 | z ( n )1 | and r r r = 1= (cid:118)(cid:117)(cid:117)(cid:116) r ( n )1 r ( n )2 | z ( n )1 | ≥ | z ( n )1 | ≥ min i | z ( n ) i | , finishing the proof. (cid:3) With the help of these two lemmas we can directly conclude Theorem 1:
Proof of Theorem 1.
For i ∈ Z consider the sequence (cid:18) r ( n ) i = | z ( n ) i − || z ( n ) i | (cid:19) n ≥ . Wehave for n ≥ i | z ( n ) i | max i | z ( n ) i | ≤ r ( n ) i ≤ max i | z ( n ) i | min i | z ( n ) i | . According to Lemmas 1.1 and 1.2, the sequence (cid:18) max i | z ( n ) i | min i | z ( n ) i | (cid:19) n ≥ is a strictly de-creasing sequence bounded from below by 1, so it converges to 1; in the sameway, the sequence (cid:18) min i | z ( n ) i | max i | z ( n ) i | (cid:19) n ≥ is a strictly increasing sequence bounded fromabove by 1, so it also converges to 1. These two results combine with (4) tolim n →∞ r ( n ) i = 1 for i = 0 , , (cid:3) Theorem 1 directly gives us the required result for the geometric element trans-formation where we assume that ∆ n is non-degenerate, that is, the vertices arepairwise disjoint: Corollary 1.1 (Elementwise convergence) . The triangle ∆ n = ( z ( n )0 , z ( n )1 , z ( n )2 ) converges for n → ∞ to an equilateral triangle. DIMITRIS VARTZIOTIS AND DORIS BOHNET
Proof of Corollary 1.2. As r ( n ) i = | z ( n ) i − z ( n ) i | converges to 1 with n → ∞ , we get thatlim | z ( n ) i | = lim | z ( n ) j | for i, j = 0 , ,
2. Therefore, the distances from the vertices tothe centroid c become equal, so that c becomes the circumcenter, and the triangleequilateral. (cid:3) Convergence of the smoothing algorithm for triangle meshes
Transformation (1) can be used to transform a mesh of triangles by combiningit at every vertex with taking the barycenter. We give the precise definition of theconsidered mesh transformation below after specifying in detail our setting.We prove in this section that any triangle mesh which does not contain too patho-logical triangles converges under the transformation to a mesh of triangles as regularas possible. Let us make precise our setting:3.1.
Preliminary notations:
Let Σ = { , . . . , N − } be a finite set of symbols.Let C = (cid:8) ∆ i = ( i , i , i ) ∈ Σ (cid:12)(cid:12) i = 0 , . . . , n − (cid:9) be a finite set of triples of symbols.We call the set C a connectivity iff for any pair ∆ i , ∆ j ∈ C there exist k ≤ n − , . . . , ∆ k ∈ C such that ∆ = ∆ i and ∆ k = ∆ j and for m = 0 , . . . , k − m and ∆ m +1 have exactly two symbols in common.Let M C : Σ → R , k (cid:55)→ x k be an injective map. We call M C a triangle meshwith connectivity C iff for any i, j = 0 , . . . , n − i (cid:54) = j the triangles defined by M C (∆ i ) := ( M C ( i ) , M C ( i ) , M C ( i )) and M C (∆ j ), counted counter clockwisely,are non-degenerate and have disjoint interior. Let denote by X C the set of meshes M C with connectivity C . Remark 1.1.
As a consequence of the definition of connectivity, the set (cid:83) n − i =0 M C (∆ i ) is arcwise connected. Definition of the mesh transformation.
We define a triangle mesh trans-formation in two steps. First, for any i = 0 , . . . , n − M C (∆ i ) = ( x i , x i , x i ) with centroid c i the triangle transformation (as definedabove) by θ i : R → R x i := ( x i , x i , x i ) (cid:55)→ ( θ i ( x i ) , θ i ( x i ) , θ i ( x i )) ,θ i ( x i ) = 23 (cid:107) x i − c i (cid:107)(cid:107) x i − c i (cid:107) ( x i − c i ) − (cid:107) x i − c i (cid:107)(cid:107) x i − c i (cid:107) ( x i − c i ) − (cid:107) x i − c i (cid:107)(cid:107) x i − c i (cid:107) ( x i − c i ) + c i θ i ( x i ) = 23 (cid:107) x i − c i (cid:107)(cid:107) x i − c i (cid:107) ( x i − c i ) − (cid:107) x i − c i (cid:107)(cid:107) x i − c i (cid:107) ( x i − c i ) − (cid:107) x i − c i (cid:107)(cid:107) x i − c i (cid:107) ( x i − c i ) + c i θ i ( x i ) = 23 (cid:107) x i − c i (cid:107)(cid:107) x i − c i (cid:107) ( x i − c i ) − (cid:107) x i − c i (cid:107)(cid:107) x i − c i (cid:107) ( x i − c i ) − (cid:107) x i − c i (cid:107)(cid:107) x i − c i (cid:107) ( x i − c i ) + c i . (5)We can now define the mesh transformation under consideration. For k = 0 , . . . , N − k denote the set of indices of the adjacent triangles at x k and we define themap Θ : X C ⊂ R N → R N x = ( x , . . . , x N − ) (cid:55)→ (Θ ( x ) , . . . , Θ N − ( x )) , x k , Θ k ( x ) ∈ R , with Θ k ( x ) = 1 | Σ k | (cid:88) m ∈ Σ k t m j ( m ) ( x k ) , (6) ONVERGENCE PROPERTIES OF A GEOMETRIC MESH SMOOTHING ALGORITHM 7 where j ( m ) ∈ , , x k inside the triangle numberedby m , so θ m j ( m ) : R → R . Remark 1.2.
The map Θ is clearly well-defined as a map from X C to R N , butnot as a map to X C : let M C (Σ) = ( x , . . . , x N − ) ∈ R N be a triangle mesh withconnectivity C . Then Θ( M C (Σ)) is not necessarily a triangle mesh M (cid:48) C . It couldhappen that the interiors of two triangles Θ(∆ i ) , Θ(∆ j ) are no longer disjoint. On the other hand note that a mesh of equilateral triangles is fixed under Θ,such that we can prove the following lemma where we call distortion of a trianglethe ratio of the shortest by the longest edge length of a triangle:
Lemma 1.3.
Well-definedness of Θ Let Θ be defined as above. Then there exists < δ < such that for any triangle mesh M C whose distortion of triangles isbounded from below by δ the image Θ( M C ) is a triangle mesh M (cid:48) C . We postpone the proof to Subsection 3.3.3.3.
Similarity group action and equivariance of Θ . Denote by Sim( R ) thefour-dimensional group of similarities of R composed by the one-dimensional group R + of scaling and the three-dimensional group Isom( R ) of isometries. This groupsnaturally acts on the set of triangles bySim( R ) × ( R ) → ( R ) ; ( g, x ) (cid:55)→ g.x = ( g.x , g.x , g.x )for any triangle x = ( x , x , x ) ∈ R where there exists A ∈ SO(2 , R ), θ ∈ R and λ ∈ R + such that g.x i = λ ( Ax i + θ ) for i = 0 , , Remark 1.3.
Two triangles x, y are similar iff there exists g ∈ Sim( R ) such that g.x = y . One can easily prove that the group Sim( R ) acts freely on the set of triangles: Lemma 1.4.
The group action as defined above is free.Proof.
Let A ∈ SO(2 , R ) , θ ∈ R and λ ∈ R + such that λ ( Ax + θ ) = x for atriangle x = ( x , x , x ) ∈ R . This implies immediately that λ = 1. So we have Ax + θ = x and Ax + θ = x . We conclude that A ( x − x ) = x − x . Thismeans that x − x is the eigenvector to an eigenvalue 1 of A . So we can concludethat A is the identity. Consequently, we get θ = 0 finishing the proof. (cid:3) The action defined above can be straightforwardly generalized to the set X C ofmeshes with connectivity C by( g, x ) ∈ Sim( R ) × ( R ) N (cid:55)→ ( g.x , g.x , . . . , g.x N − )for any mesh M C (Σ) =: x ∈ ( R ) N . This action is certainly free as well. Remark 1.4.
One could think of defining the similarity group action on a meshseparately on every triangle. But in fact, the connectivity as defined above forcesthat the same group element acts simultaneously on each triangle of the mesh. Sothe group action defined above is the only one in accordance with the given definitionof a mesh.
As a consequence we can list the following properties of the group action:(1) Every Sim( R )-orbit is a four-dimensional smooth submanifold in R N .(2) One computes immediately that the mesh transformation Θ is equivariant under the group action of Sim( R ), that is(a) For every M C inside the domain of Θ the image Θ( M C ) lies as well inthe domain. DIMITRIS VARTZIOTIS AND DORIS BOHNET (b)Θ( g.M C ) = g. Θ( M C ) for any M C ∈ X C and g ∈ Sim( R ) . (see e.g. [4] where important properties for equivariant dynamical systemsare proved).(3) For any h ∈ Sim( R ) one computes for the Jacobian matrix of Θ for any M C ∈ X C Θ = h − ◦ Θ ◦ h ⇐ D Θ M C = Dh − ◦ D Θ h.M C ◦ Dh, and as h, h − are linear maps one gets D Θ M C = h − ◦ D Θ h.M C ◦ h. Thanks to these properties we can reduce the question of global convergence to thefollowing: Let M C be a fixed point of Θ, that is Θ( M C ) = M C , then the wholegroup orbit Λ := Sim( R ) .M C is fixed and Λ is consequently a four-dimensionalsubmanifold of fixed points.Now we can prove Lemma 1.3: Proof of Lemma 1.3.
Let M eq be an equilateral mesh, then we have Θ( M eq ) = M eq .On the other hand, the domain of Θ is clearly an open subset of X C . By thecontinuity of Θ, there exists an open set U of triangles sufficiently close to M eq suchthat Θ is well-defined on U . The equivariance of Θ implies that Θ is well-definedon the group orbit Sim( R ) . U of U which contains all meshes whose distortion oftriangles is bounded by some 0 < δ < (cid:3) To study the convergence in a neighborhood of the fixed point M C it is enoughto study the dynamics of Θ in a neighborhood of Λ thanks to the following lemma: Lemma 1.5.
Let x = M C ∈ R N be a fixed point of Θ and Λ := Sim( R ) .M C itsgroup orbit. If there exists a D Θ -invariant decomposition of the tangent bundle at Λ T R N | Λ = T Λ ⊕ E s such that (cid:107) D Θ | E s (cid:107) < , then there exists a unique family F s of injectively C r -immersed submanifolds F s ( x ) such that x ∈ F s ( x ) and F s ( x ) is tangent to E sx atevery x ∈ Λ . This family is Θ -invariant, that is, Θ( F s ( x )) = F s (Θ( x )) , and themanifolds F s ( x ) are uniformly contracted by some iterate of Θ .That family actually forms a foliation of a neighborhood of Λ . This Lemma is an immediate application of the invariant manifold theorem byHirsch,Pugh and Shub, cited and proved for example in [1, Th.B7, p.293]. Thespectrum spec( D Θ | Λ ) tangent to the group orbit contains four eigenvalues equal to1. Consequently, we have the following direct corollary of Lemma 1.5: Corollary 1.2.
Let M C ∈ R N be a fixed point of Θ and Λ := Sim( R ) .M C itsgroup orbit. If every eigenvalue of the Jacobian matrix D Θ which is not containedin spec( D Θ | Λ ) has an absolute value strictly smaller than , then Λ is an attractorand Θ n ( M ) converges uniformly at exponential rate to one point in Λ for n → ∞ and for any triangle mesh M sufficiently close to Λ . So as a consequence of this corollary, it is enough to study the spectrum of D Θat a fixed point and to prove that the absolute value of all eigenvalues except fromfour is strictly smaller than one. Nevertheless, this is still a difficult task as it willbecome obvious in the following. We start with the easiest cases gaining more andmore complexity:
ONVERGENCE PROPERTIES OF A GEOMETRIC MESH SMOOTHING ALGORITHM 9
Convergence for particular cases.
Case 1: A single triangle.
We start with the easiest case of a mesh which con-sists of a single triangle, so in fact, we study the global convergence of the previouslydefined triangle transformation θ : R → R on a triangle x = ( x , x , x ) ∈ ( R ) in more details and using the new setting above. Let x eq ∈ R be an equilateraltriangle, then θ ( x eq ) = x eq . The equivariance of θ under the group of similaritiesprovokes that θ (Λ eq ) = Λ eq , where Λ eq = (cid:8) x ∈ R (cid:12)(cid:12) g ∈ Sim( R ) : x = g.x eq , (cid:9) . is a 4-dimensional submanifold of R which is θ -invariant, that is, θ (Λ eq ) ⊂ Λ eq .Following Corollary 1.2 we compute the derivative Dθ x eq of θ at x eq ∈ Λ eq . TheJacobian matrix is the same matrix J for any x eq ∈ Λ eq : Dt ( x eq ) = A B CC A BB C A =: J where A = (cid:32) − √ √ (cid:33) , B = (cid:32) − √ √ (cid:33) , C = (cid:32) √ − √ (cid:33) . (7)Remark that J is a circulant block matrix. Further, J is conjugate to the blockdiagonal matrix ( R π/ , id R ) where R π/ : R → R is a rotation by π/
3. Accord-ingly, there exists a 2-dimensional subspace E s spanned by the eigenvectors v , v corresponding to the two eigenvalues (cid:54) = 1. There exists c > v ∈ E sx , x ∈ Λ eq , one has (cid:107) Dt x v (cid:107) ≤ c (cid:107) v (cid:107) . The four eigenvectors v , . . . , v corresponding to eigenvalues 1 span the tangentspace of Λ eq . So – applying Corollary 1.2 – the invariant set Λ eq is an attractor for θ . Hence, there exists a neighborhood U eq ⊃ Λ eq such that every x ∈ U eq convergesto Λ eq under iterates of θ , that is,dist( θ n x, Λ eq ) → , n → ∞ . Taking into account Theorem 1 one concludes that Λ eq is a global attractor . Remark 1.5.
If one considers the orbit space of the free group action
Isom( R ) on R by identifying similar triangles, one observes that this space is the two-dimensional projective space P ( R ) . By the observation above the triangle trans-formation passes to a well defined map on this quotient space: R p (cid:15) (cid:15) t (cid:47) (cid:47) R p (cid:15) (cid:15) P ( R ) t (cid:47) (cid:47) P ( R ) The attractor Λ eq projects to a globally attracting fixed point on P ( R ) . By theequivariance of the transformation θ , the two-dimensional stable set tangent to E s passes also to a well-defined two-dimensional set in the orbit space reflecting theattraction of the fixed point. Case 2: mesh of six equilateral triangles.
Let Σ = { , . . . , } and C = { (0 , , , (0 , , , (0 , , , (0 , , , (0 , , , (0 , , } be the connectivity and denoteby x eq = ( x , . . . , x ) ∈ ( R ) the mesh of six equilateral triangles. The grouporbit Λ eq under the similarity group action is – exactly as above – a 4-dimensionalsmooth submanifold of R , and every mesh x eq ∈ Λ eq is certainly a fixed point ofthe mesh transformation Θ. We compute – with the notations above – the Jacobianmatrix of Θ for x eq ∈ Λ eq as D Θ( x eq ) = A ( B + C ) T ( B + C ) T ( B + C ) T ( B + C ) T ( B + C ) T ( B + C ) T B A B T C T B C T A B T B C T A B T B C T A B T B B T C T A B T B C T A . (8)The matrix D Θ( x eq ) has four eigenvalues 1 whose eigenvectors span the 4-dimensionaltangent space of Λ eq . Further, we have five pair of complex conjugate eigenvalues λ , λ , . . . , λ with absolute values | λ i | ∈ [0 . , . x eq ∈ Λ eq splits into a ten dimensional space E s ( x eq ) spanned bythe eigenvectors v , v , . . . , v and a 4-dimensional eigenspace T x eq Λ eq of the equiv-alence relation: T Λ eq R = E s (Λ eq ) ⊕ T Λ x eq . So we can apply Corollary 1.2 and conclude that Λ eq is a local attractor, andconsequently, there exists a neighborhood Λ eq ⊂ U eq ⊂ R such that every mesh x = M C ∈ U eq converges uniformly to one mesh x eq under Θ:dist(Θ n ( x ) , x eq ) → n →∞ x ∈ U eq . Remark 1.6.
In contrast to the case of the triangle transformation we cannot provethat Λ eq is a global attractor: one observes numerically that D Θ( x ) for x ∈ X C might have eigenvalues of absolute value > , that is, there are directions in which x is expanded. Numerical tests show, that after one or two iterations of Θ , x comessufficiently close to x eq such that it converges uniformly to x eq . Case 3: simple meshes. let Σ = { , . . . , N − } be a set of N symbols. Wecall a connectivity C N - simple iff all triples ( i , i , i ) ∈ C has a common symbol.We call M C a N -simple mesh iff its connectivity C is N -simple.Above, we consider – in this terminology – a 6-simple mesh. For a N -simple mesh,one fixed point of Θ is the mesh defined by the vertices x = (0 ,
0) and x k − =(cos(2 kπ/ ( N − , sin(2 kπ/ ( N − k = 2 , . . . , N . Let denote the similaritygroup orbit of this mesh by Λ eq,N . We can then numerically compute the Jacobianmatrix for x eq and show their spectra in Figure 2 for N = 4 , . . . , N,eq is for 4 < N ≤
11 a local attractor. Further,Figure 2 seems to suggest that for N ∈ [5 , eq,N is attractingin a quite large region.3.4.4. Case 4: mesh of equilateral triangles.
Let C be a connectivity such that everyinner vertex has exactly six neighboring vertices, and consider the previously definedset X C of meshes with this connectivity. Denote by N i the indices of inner verticesand by N b the indices of boundary vertices. Let x eq = ( x , . . . , x N − ) ∈ ( R ) N bethe mesh of equilateral triangles. Exactly as above, we consider the whole group ONVERGENCE PROPERTIES OF A GEOMETRIC MESH SMOOTHING ALGORITHM 11
102 4 6 8 123 5 7 9 11010.20.40.60.81.21.41.6 N−simple mesh ab s o l u t e v a l ue o f e i gen v a l ue s f o r J a c ob i an m a t r i x
102 4 6 8 123 5 7 9 11010.20.40.60.81.21.41.6 N−simple mesh ab s o l u t e v a l ue o f e i gen v a l ue s f o r J a c ob i an m a t r i x Figure 2.
Plot of the spectrum of DT N ( x ) where T N is the meshtransformation of a N -simple meshes and x runs through 50 ran-domly generated N -simple meshes. On the right, the spectrum forthe regular N -simple mesh is depicted. Note that the equilateral3-simple mesh is not an attracting point, but a saddle point .orbit Λ eq ⊂ X C of the similarity group. Then we compute the Jacobian matrix D Θ( x eq ) of the mesh transformation (6) at x eq ∈ Λ eq : D Θ( x eq ) = (cid:18) ∂ Θ k ∂x l (cid:19) k,l =0 ,...,N − , ∂ Θ k ∂x l ∈ R × ∂ Θ k ∂x l = A if k = l∂ Θ k ∂x l = 16 ( B + C ) T if l ∈ Σ k , k ∈ N i ,∂ Θ k ∂x l = 12 C if l ∈ Σ k , l, k ∈ N b ,∂ Θ k ∂x l = 12 B if l ∈ Σ k , l ∈ N i , k ∈ N b ,∂ Θ k ∂x l = 0 if l / ∈ Σ k . (9)After various computations on different equilateral meshes we conjecture the fol-lowing: Conjecture 1.1.
For any equilateral mesh x the Jacobian matrix of Θ at x haseigenvalues of absolute value < except from exactly four. In particular, the grouporbit Sim( R ) .x of the mesh x is an attractor. Further generalization.
One could again study the jacobian matrix of Θ atany fixed point x . But things get much more complicated, because the matricescould not be expressed in a simple way. We conjecture the following, where ˜ X C isthe quotient space X C / Sim( R ) of the group action of Sim( R ): Conjecture 1.2.
For any ≤ N < ∞ and any connectivity C with cardinality N the following is true: there exists a metric (cid:107) (cid:107) X on the quotient space ˜ X C such thatthe map ˜Θ induced on ˜ X C is strictly contracting on its domain with respect to thismetric, that is (cid:107) T ( M C ) − T ( M (cid:48) C ) (cid:107) X ≤ λ (cid:107)M C − M (cid:48) C (cid:107) X for any two meshes M C , M (cid:48) C ∈ ˜ X C . This would imply in particular that any fixed point ˜ x ∈ ˜ X C is an attractor. Remark 1.7. (1)
In Figure 3 we show the absolute value of the six eigenvaluesfor 700 randomly generated triangles. This figure stresses also the factthat for not too distorted triangles the triangle transformation is strictlycontracting transverse to the normally hyperbolic invariant set characterizedby the four eigenvalues equal to one. (2)
In Figure 5 we computed the norm of the Jacobian of the mesh transforma-tion of a mesh of triangles (shown in the left picture) in relation to thematrix norm of the Jacobian for the most regular mesh of triangles. Oneobserves in the right picture how the matrix norm approaches the optimalmatrix norm as the quality of the triangle mesh approaches its optimum. ab s o l u t e v a l ue o f e i gen v a l ue lambda1lambda2lambda3lambda4lambda5lambda6 Figure 3.
Plot ofthe absolute valueof eigenvalues independence of thetriangle quality of 700randomly generatedtriangles.
102 4 6 81 3 5 7 900.20.40.10.30.050.150.250.350.45 iteration of triangle transformation de v i a t i on o f F r oben i u s no r m Figure 4.
Plotof the deviation | (cid:13)(cid:13) dt i (cid:13)(cid:13) − (cid:107) dt e (cid:107) | of theFrobenius norm of t i , i = 1 , . . . ,
10 fromthe Frobenius norm (cid:107) dt e (cid:107) correspondingto an equilateraltriangle.Outlook: The proof should be easily generalized for triangle meshes defined on Rie-mannian surfaces, that is, – with the notations above – the triangle mesh is definedby M C : Σ → S, k (cid:55)→ x k ∈ S where S is a Riemannian surface such that everytriangle M (∆) lies inside one chart neighborhood.The techniques developed in this proof could also be adaptable to similar geometricmesh transformations.In [13], we model the triangle transformation above by system of linear differentialequations which could be seen as the description of coupled damped oszillations.This model provides another explanation why the transformation converges to aequilateral triangle. One could think of the mesh transformation as the discretiza-tion of the solution of a system of coupled damped oszillations which are driven byeach other antagonizing the damping.4. Short discussion of implementation and numerical results
We do not focus in this article on the application of our algorithm, so the follow-ing discussion is kept very brief and should be treated as a motivation to explorefurther the practical possibilities of the presented algorithm in the future. We haveimplemented the method as it is described above in Section 3 inside the open sourcesoftware Scilab 5.4.1. The method could be equally well directly implemented in
ONVERGENCE PROPERTIES OF A GEOMETRIC MESH SMOOTHING ALGORITHM 13
Figure 5.
Plot of (cid:107) dT (cid:107) and mean mesh quality during iterationsfor a mesh of 7 triangles together with the same values for anoptimal mesh of 7 triangles with inner angle 2 π/ C . For an industrial usage this is strongly preferable to make it more efficient.We tested the method for a randomly generated triangulation of the unit square.See below in Figure 6 how the mesh converges to a mesh of quite equilateral trian-gles in very few iterations. As quality measure q ∆ we used the ratio of minimal to
01 0.20.40.60.8 0.10.30.50.70.9 010.20.40.60.80.10.30.50.70.9X YZ 01 0.20.40.60.8 0.10.30.50.70.9 010.20.40.60.80.10.30.50.70.9X YZ
Figure 6.
A randomly generated triangulation of the unit squareat the beginning and after 10 iterations colored depending on theirquality measure q ∆ ∈ (0 , q ∆ = min i,j =1 (cid:107) x i − x j (cid:107) max i,j =1 (cid:107) x i − x j (cid:107) . The quality measurefor a triangle mesh V = (∆ , . . . , ∆ | V |− ) is then the mean of the quality measure q ∆ for every triangle ∆ ∈ V : q V = | V | (cid:80) ∆ ∈ V q ∆ . The mesh we smoothed in Figure 6 consists of 450 triangle elements. In Figure 8bbelow we show how the number of elements with a certain quality measure developsover iterating the mesh and how the mean quality improves.The smoothing algorithm works equally well for tetrahedra by applying the smooth-ing algorithm to the triangular faces. We display in Figure 7 the cube [0 , cut at x = 0 . T = ( x , x , x , x ) with x i ∈ R be a tetrahedron. As quality measure q T for a tetrahedron T we use the mean ratio quality measure which is defined as (a) Initial mesh: q V = 0 . (b) q V = 0 . Figure 7.
Application to a tetrahedral mesh of 5318 elements ofthe unit cube (for x < . q T ( T ) = 3 det( S ) / trace( S t S ) , S = D ( T ) W, where D ( T ) = ( x − x , x − x , x − x ) , W = / / √ / √ /
60 0 (cid:112) / . As quality measure for a tetrahedral mesh V = ( T , . . . , T | V |− ) we used the meanquality measure of every element: q V = | V | (cid:80) T ∈ V q T ( T ) . In Figure 8a, one canobserve how the quality measure of the mesh of Figure 7 improves. (a)
Improvement of the cube meshelement quality for the mesh in Fig-ure 7 (b)
Square quality element measure q ∆ before and after the smoothingof the mesh in Figure 6. Concluding remarks
Generalization to polygonal meshes.
Any polygon can be transformed inthe exactly analogous way as the triangle above. Let P = ( x (0)0 , . . . , x (0) k − ) be aconvex k -gon with x (0) i ∈ R with its centroid in the origin. Then we can define atransformation in the following way recursively: x ( n +1) i = k − k r ( n ) i x ( n ) i − k k − (cid:88) j =0 ,j (cid:54) = i r ( n ) j x ( n ) j , r ( n ) i = (cid:13)(cid:13)(cid:13) x ( n ) i − (cid:13)(cid:13)(cid:13) (cid:13)(cid:13)(cid:13) x ( n ) i (cid:13)(cid:13)(cid:13) . Remark that the centroid is kept in the origin througout the transformation. Butthe iterated polygon P ( n ) = ( x ( n )0 , . . . , x ( n ) k − ) does not necessarily converge for n →∞ to a polygon with equal distances (cid:107) x i (cid:107) = (cid:107) x j (cid:107) , i, j = 0 , . . . , k −
1, i.e. its
ONVERGENCE PROPERTIES OF A GEOMETRIC MESH SMOOTHING ALGORITHM 15 centroid coincides with its circumcenter. Consider for example a quadrilateral Q (0) with (cid:107) x (cid:107) = (cid:107) x (cid:107) and (cid:107) x (cid:107) = (cid:107) x (cid:107) . Then Q (0) = Q (2) is two-periodic but do notconverge. So the transformation has not a globally attracting fixed point for allinitial polygons. Also observe the following: while a triangle is regular if and onlyif the distances of its vertices to its centroid is equal, this is not the case for otherpolygons where it is just a necessary, but not a sufficient condition.Accordingly, the transformation cannot be directly used for a smoothing algo-rithm for polygonal meshes without further adaption.But nevertheless, the transformation can be used to smooth any polygonal meshby subdividing every polygon into triangles and then applying the transformationto every triangle.5.2. Outlook.
One easily detects the following shortcomings of the presented smooth-ing method which are open for future research:
Global convergence:
We only prove the global convergence for a compactsubset of triangle meshes which exclude triangles close to degenerate ones.By changing the transformation a bit – with regard to the estimates wederive during the proof – such that the transformation is integrable, thatis, the gradient of a function, and consequently the jacobian matrix normal,one could obtain better bounds and therefore extend the convergence resultto a greater subset of meshes.
Performance:
It was not the primary objective of this article to provide anefficient implementation, but to analyze the underlying mathematics. Sowe have to admit that each iteration step is numerically quite long in thepresent implementation. But if directly implemented inside C , we shouldattain comparable run times as for GETMe. As the algorithm relocatesseparately every vertex, it is open for an application of parallel computingtechniques. Polygonal/ polyhedral meshes:
Using the duality of certain polygons/polyhedra to each other we hope to be able to adapt the current trans-formation to quadrilateral and hexahedral meshes. This is a current topicof our research.
References [1] Christian Bonatti, Lorenzo Diaz, and Marcelo Viana.
Dynamics beyond uniform hyperbolicity:a global geometric and probabilistic perspective . Springer, 2005.[2] Graham Carey.
Computational grids: generation, adaptation and solution strategies . Taylorand Francis, 1998.[3] David A. Field. Laplacian smoothing and Delaunay triangulations.
Communications in Ap-plied Numerical Methods , 4(6):709–712, 1988.[4] Michael Field. Equivariant dynamical systems.
Trans. Amer. Math. Soc. , 259(1):185–205,1980.[5] Lori A. Freitag. On combining Laplacian and optimization-based mesh smoothing techniques.In
Trends in Unstructured Mesh Generation , pages 37–43, 1997.[6] Lori A. Freitag and Patrick M. Knupp. Tetrahedral element shape optimization via the ja-cobian determinant and condition number. In
Proceedings of the 8th International MeshingRoundtable , pages 247–258. Sandia National Laboratory, 1999.[7] Pascal J. Frey and Paul-Louis George.
Mesh Generation . Hermes Science Publishing, 2000.[8] Leonard Herrmann. Laplacian-isoparametric grid generation scheme.
Journal of the Engi-neering Mechanics Division , 102(5):749–756, 1976.[9] Mike Holder and Charles Karr. Quadrilateral mesh smoothing using a steady state geneticalgorithm. In
Genetic and Evolutionary Computation GECCO 2003 , volume 2724 of
LectureNotes in Computer Science , pages 2400–2401. Springer Berlin Heidelberg, 2003.[10] Patrick M. Knupp. Algebraic mesh quality metrics.
SIAM Journal on Scientific Computing ,23(1):193–218, 2001. [11] Dimitris Vartziotis. General transformations – regularization and symmetry. unpublishedmanuscript, 2013.[12] Dimitris Vartziotis, Theodoros Athanasiadis, Iraklis Goudas, and Joachim Wipper. Meshsmoothing using the geometric element transformation method.
Comput. Methods Appl.Mech. Engrg. , 197(45-48):3760–3767, 2008.[13] Dimitris Vartziotis and Doris Bohnet. A geometric triangle and tetrahedral mesh smoothingalgorithm related to damped oszillations. Preprint, 2014.[14] Dimitris Vartziotis and Benjamin Himpel. Efficient and global optimization-based smoothingmethods for mixed-volume meshes. In Josep Sarrate and Matthew Staten, editors,
Proceed-ings of the 22nd International Meshing Roundtable , pages 293–311. Springer InternationalPublishing, 2014.[15] Dimitris Vartziotis, Joachim Wipper, and Bernd Schwald. The geometric element transforma-tion method for tetrahedral mesh smoothing.
Comput. Methods Appl. Mech. Engrg. , 199(1-4):169–182, 2009.[16] A. Egemen Yilmaz and Mustafa Kuzuoglu. A particle swarm optimization approach for hex-ahedral mesh smoothing.
Int. J. Numer. Meth. Fluids , 60(1):55–78, 2009.
TWT GmbH Science & Innovation, Department for Mathematical Research & Ser-vices, Ernsthaldenstr. 17, 70565 Stuttgart, Germany;NIKI Ltd. Digital Engineering, Research Center, 205 Ethnikis Antistasis Street, 45500Katsika, Ioannina, Greece
E-mail address : [email protected] TWT GmbH Science & Innovation, Department for Mathematical Research & Ser-vices, Ernsthaldenstr. 17, 70565 Stuttgart, Germany
E-mail address ::