A perimeter-decreasing and area-conserving algorithm for surface diffusion flow of curves
aa r X i v : . [ m a t h . NA ] J a n A perimeter-decreasing and area-conserving algorithmfor surface diffusion flow of curves
Wei Jiang a , Buyang Li b, ∗ a School of Mathematics and Statistics & Hubei Key Laboratory of Computational Science,Wuhan University, Wuhan 430072, P. R. China. b Department of Applied Mathematics, The Hong Kong Polytechnic University, Hong Kong.
Abstract
A fully discrete finite element method, based on a new weak formulation and anew time-stepping scheme, is proposed for the surface diffusion flow of closedcurves in the two-dimensional plane. It is proved that the proposed method canpreserve two geometric structures simultaneously at the discrete level, i.e., theperimeter of the curve decreases in time while the area enclosed by the curve isconserved. Numerical examples are provided to demonstrate the convergence ofthe proposed method and the effectiveness of the method in preserving the twogeometric structures.
Keywords: surface diffusion flow, area conservation, perimeter decrease,parametric, weak formulation, time stepping, finite element method
1. Introduction
This article concerns the numerical approximation to the surface diffusionflow of closed curves in the two-dimensional plane, i.e., the evolution of a curveΓ[ X ( · , t )] = { X ( ξ, t ) : ξ ∈ I } , t ∈ [0 , T ] , ✩ This work is supported by the National Natural Science Foundation of China (project no.11871384) and the Hong Kong Research Grants Council (GRF project no.: 15300920). ∗ Corresponding author
Email addresses: [email protected] (Wei Jiang), [email protected] (Buyang Li)
Preprint submitted to Journal of Computational Physics February 2, 2021 etermined by a parametrization X ( · , t ) : I → R , where I = [0 ,
1] is theperiodic unit interval (one-dimensional torus which identifies 0 and 1), satisfyingthe fourth-order geometric evolution equation ∂ t X = − [ ∂ ( κ · ν )] ν , κ = ∂ Γ τ , on I × [0 , T ] , (1)where κ = κ ( ξ, t ), ν = ν ( ξ, t ) and τ = τ ( ξ, t ) are the curvature vector, inwardunit normal vector and unit tangential vector of the curve Γ[ X ( · , t )] at the point X ( ξ, t ). The notation ∂ Γ represents tangential differentiation along the curve Γ,defined by ∂ Γ v ( ξ, t ) := ∂ ξ v ( ξ, t ) | ∂ ξ X ( ξ, t ) | for a function v defined on I × [0 , T ] . The parametric equation (1) for surface diffusion flow was proposed byMullins in 1957 for modeling the evolution of microstructure in polycrystallinematerials [27]. It was generalized to include anisotropic effects of crystallinefilms in [12, 20, 26, 30]. These models play important roles in various applica-tions in materials science and solid-state dewetting; see [3, 21–25, 29, 31] andthe references therein.It is well-known that the surface diffusion flow described by (1) is the H − gradient flow of the perimeter functional (see [12, 28]) | Γ[ X ( · , t )] | = Z I | ∂ ξ X ( ξ, t ) | d ξ. As a result, the surface diffusion flow has two geometric structures:(i) The perimeter of the evolving curve decreases in time;(ii) The area enclosed by the curve is conserved in the evolution.These geometric structures are also desired in the numerical approximation ofsurface diffusion flow.Many numerical methods have been developed for simulating the evolutionof a curve/surface in 2D/3D governed by surface diffusion flow and related2eometric evolution equations. For the non-parametric equations describing ax-ially symmetric surface diffusion flow, Coleman, Falk & Moakher [13] proposedspace-time finite element methods (FEMs) to preserve the perimeter decreaseand area conservation properties, and Deckelnick, Dziuk & Elliott [14] provedconvergence of semidiscrete FEMs. For the non-parametric equations describ-ing surface diffusion flow of graphs, variational formulation, error analysis andnumerical simulation were done by B¨ansch, Morin & Nochetto [1] and Deckel-nick, Dziuk & Elliott [15] for semidiscrete and fully discrete FEMs, respectively.In the axially symmetric and graph cases, both semi-discrete and fully discretemethods can preserve the perimeter decrease and area conservation.For general parametric surface diffusion flow, the underlying numerical meth-ods are the parametric FEMs developed by Dziuk [17, 18], Deckelnick, Dziuk &Elliott [16], B¨ansch, Morin & Nochetto [2], Dziuk & Elliott [19], Barrett, Garcke& N¨urnberg [5, 7], etc. In particular, B¨ansch, Morin & Nochetto [2] proposedthe first parametric FEM for the surface diffusion flow. The method preservesthe decrease of perimeter (see [2, Theorem 2.1]), while the area enclosed bythe curve may change in time. Parametric FEMs based on novel variationalformulations were developed in [5, 7, 9, 11]. These methods yield good meshdistribution and unconditional stability. The generalization to anisotropic sur-face diffusion flows and axisymmetric geometry cases was made in [6, 8, 10].These methods can preserve the decrease of perimeter in 2D (or decrease ofarea in 3D) in the fully-discrete finite element scheme, and the area conserva-tion in 2D (or volume conservation in 3D) only in the semi-discrete FEM. Morerecently, Jiang et al. developed parametric FEMs for simulating solid-statedewetting problems described by surface diffusion flow and contact line migra-tion in 2D [4] and 3D [32], without considering the perimeter/area decrease andarea/volume conservation properties in the discrete level.Overall, existing numerical methods for parametric surface diffusion flow of-ten yield curves/surfaces with decreasing perimeter/area, while the area/volumeenclosed by the curve/surface is conserved only in the semi-discrete FEMs. Inthis paper, we construct a fully discrete parametric FEM for surface diffusion3ow of closed curves, based on a new weak formulation and a new time-steppingscheme, to preserve the two geometric structures simultaneously.The rest of paper is organized as follows. In section 2, we propose a newweak formulation for solving problem (1), then a parametric FEM is used todiscretize the weak formulation and results in a nonlinear system of equations.In section 3, we rigorously prove that the numerical scheme is area-preservingand perimeter-decreaseing. In section 4, we present some numerical simulationsto demonstrate the accuracy and efficiency of the proposed method. Finally, wedraw some conclusions in section 5.
2. The numerical method
In this section, we introduce a new variational formulation of the surfacediffusion flow equation (1) and then propose fully discrete FEMs correspondingto the variational forms.
Let 0 = ξ < ξ · · · < ξ M = 1 be a quasi-uniform partition of the interval I (where we identify ξ and ξ M ), which is divided into M subintervals labeled as I j := [ ξ j − , ξ j ], j = 1 , · · · , M . We define the following finite element spaces ofpiecewise linear functions subject to the partition, i.e., V h = (cid:8) ϕ ∈ C ( I ; R ) : ϕ (0) = ϕ (1) , ϕ | I j ∈ P , j = 1 , , · · · , M (cid:9) , (2) V h = (cid:8) ψ ∈ C ( I ; R ) : ψ (0) = ψ (1) , ψ | I j ∈ P × P , j = 1 , , · · · , M (cid:9) , (3)where P denotes the space of polynomials of degree ≤ t n = nτ , n = 0 , , , . . . , N , be a uniform partition of the time interval[0 , T ] with stepsize τ = T /N . For each n = 1 , . . . , N , we denote by Γ nh apiecewise linear curve determined by a sequence of non-coincident nodes x nj ∈ R , j = 1 , · · · , M in the counter clockwise order, approximating the closed curveΓ[ X ( · , t n )]. Correspondingly, the nodal vector x n = ( x n , · · · , x nM ) uniquelydetermines a finite element function X nh ∈ V h which parametrizes Γ nh , satisfying X nh ( ξ j ) = x nj , j = 1 , · · · , M.
4e assume that X n − h ∈ V h is given and look for X nh ∈ V h by an implicittime-stepping scheme defined below.In the time interval [ t n − , t n ], we define the intermediate curvesΓ h ( t ) = { X h ( ξ, t ) : ξ ∈ [0 , } with parametrization X h ( ξ, t ) = t n − tτ X n − h ( ξ ) + t − t n − τ X nh ( ξ ) , t ∈ [ t n − , t n ] . (4)Then for a fixed ξ ∈ I the material point X ( ξ, t ) on the intermediate curve Γ h ( t )moves with velocity v nh ( ξ ) = X nh ( ξ ) − X n − h ( ξ ) τ , (5)which is a finite element function in V h independent of time.The unit normal and tangential vectors at the X h ( ξ, t ) ∈ Γ h ( t ) are given by ν h ( ξ, t ) = ∂ ξ X h ( ξ, t ) ⊥ | ∂ ξ X h ( ξ, t ) | and τ h ( ξ, t ) = ∂ ξ X h ( ξ, t ) | ∂ ξ X h ( ξ, t ) | , (6)respectively, where ⊥ denotes rotation of a vector by an angle π/ We decompose the curvature vector into its normal and tangential compo-nents, separately, i.e., κ = p ν + q τ , where p and q are scalar functions (in fact, q = 0 for the continuous problem).With these notations, we rewrite equation (1) as ∂ t X · ν = − ∂ p, (7a) ∂ t X · τ = 0 , (7b) p ν + q τ = ∂ Γ τ . (7c)5et H ( I ) = { φ ∈ H (0 ,
1) : φ (0) = φ (1) } and consider the following weakformulation of (7): find X ( · , t ) ∈ H ( I ) × H ( I ) and p ( · , t ) , q ( · , t ) ∈ H ( I ) suchthat the equations Z I ∂ t X ( · , t ) · ( φ ν ) dΓ( t ) = Z I ∂ Γ p ∂ Γ φ dΓ( t ) , ∀ φ ∈ H ( I ) , (8a) Z I ∂ t X ( · , t ) · ( ϕ τ ) dΓ( t ) = 0 , ∀ ϕ ∈ H ( I ) , (8b) Z I ( p ν + q τ ) · ψ dΓ( t ) = Z I ∂ Γ τ · ψ dΓ( t ) , ∀ ψ ∈ H ( I ) × H ( I ) , (8c)hold for all t ∈ (0 , T ] under an initial condition X ( ξ,
0) = X ( ξ ), where X is agiven parametrization of the curve at time t = 0, anddΓ( t ) := | ∂ ξ X ( ξ, t ) | d ξ. On the intermediate curve Γ h ( t ), we approximate the curvature at the point X h ( ξ, t ) by a function κ h ( ξ, t ) = p h ( ξ ) ν h ( ξ, t ) + q h ( ξ ) τ h ( ξ, t ) , (9)with some finite element functions p h , q h ∈ V h independent of time when t ∈ [ t n − , t n ].By using the discrete velocity and curvature defined by (5) and (9), we lookfor X nh ∈ V h and p nh , q nh ∈ V h satisfying the following equations: Z t n t n − Z I v nh · ( φ h ν h ) dΓ h ( t )d t = Z t n t n − Z I ∂ Γ p nh ∂ Γ φ h dΓ h ( t )d t, ∀ φ h ∈ V h , (10a) Z t n t n − Z I v nh · ( ϕ h τ h ) dΓ h ( t )d t = 0 , ∀ ϕ h ∈ V h , (10b) Z t n t n − Z I ( p nh ν h + q nh τ h ) · ψ h dΓ h ( t )d t = M X j =1 ψ h ( ξ j ) · Z t n t n − [ τ h ] j d t, ∀ ψ h ∈ V h , (10c)where dΓ h ( t ) := | ∂ ξ X h ( ξ, t ) | d ξ and[ τ h ] j = lim ε → ( τ h ( ξ j + ε, t ) − τ h ( ξ j − ε, t ))6enotes the jump of tangential vector at the node X h ( ξ j , t ). The right side of(10c) is obtained by using the identity ∂ ξ τ h ( ξ, t ) = M X j =1 [ τ h ] j δ ( ξ − ξ j )for the piecewise constant function τ h , where δ ( ξ − ξ j ) denotes the Dirac deltafunction centered at the point ξ j .If we denote ν nh ( ξ ) = Z t n t n − ν h ( ξ, t ) | ∂ ξ X h ( ξ, t ) | d t = Z t n t n − ( ∂ ξ X h ( ξ, t )) ⊥ d t, (11) τ nh ( ξ ) = Z t n t n − τ h ( ξ, t ) | ∂ ξ X h ( ξ, t ) | d t = Z t n t n − ∂ ξ X h ( ξ, t ) d t, (12) ρ nh ( ξ ) = Z t n t n − | ∂ ξ X h ( ξ, t ) | − d t, (13)then (10) can be equivalently as: find X nh ∈ V h and p nh , q nh ∈ V h such that Z I X nh − X n − h τ · ( φ h ν nh ) d ξ = Z I ∂ ξ p nh ∂ ξ φ h ρ nh d ξ, ∀ φ h ∈ V h , (14a) Z I X nh − X n − h τ · ( ϕ h τ nh ) d ξ = 0 , ∀ ϕ h ∈ V h , (14b) Z I ( p nh ν nh + q nh τ nh ) · ψ h d ξ = M X j =1 ψ h ( ξ j ) · Z t n t n − [ τ h ] j d t, ∀ ψ h ∈ V h . (14c)The numerical scheme (14) together with (11)-(13) is a nonlinear system ofequations, which can be solved by using the following Newton’s iteration. In this subsection, we prove that the proposed numerical scheme (14), orits equivalent form (10), can preserve the two geometric structures (i) and (ii)mentioned in the introduction section.Let L ( t n ) be the perimeter of the piecewise linear curve Γ h ( t n ), and let A ( t n )be the area of the polygonal region enclosed by Γ h ( t n ). Theorem 1.
The solution of (14) has the following properties: (1) A ( t n ) = A ( t n − ) , L ( t n ) ≤ L ( t n − ) .Proof. We use the equivalent formulation (10) to prove these properties. Byintegrating the identity dd t A ( t ) = − Z I v nh · ν h dΓ h ( t ) , in time for t ∈ [ t n − , t n ] and setting φ h = 1 in (10a), we obtain A ( t n ) − A ( t n − ) = − Z t n t n − Z I v nh · ν h dΓ h ( t )d t = − Z t n t n − Z I ∂ Γ p nh · ∂ Γ h ( t )d t = 0 . This proves that the area enclosed by the curve is conserved during the evolution.We use the following expression for the perimeter of the piecewise linearcurve Γ h ( t ): L ( t ) = M X j =1 Z ξ j ξ j − | ∂ ξ X h ( ξ, t ) | d ξ. Differentiating this expression in time yieldsdd t L ( t ) = M X j =1 Z ξ j ξ j − ∂ t | ∂ ξ X h ( ξ, t ) | d ξ = M X j =1 Z ξ j ξ j − ∂ ξ X h ( ξ, t ) | ∂ ξ X h ( ξ, t ) | · ∂ ξ ∂ t X h ( ξ, t )d ξ = M X j =1 Z ξ j ξ j − τ h · ∂ ξ v h d ξ = − M X j =1 [ τ h ] j · v h ( ξ j ) , where we have used integration by parts on each subinterval [ ξ j − , ξ j ]. Byintegrating the equality above with respect to time from t n − to t n , we obtain L ( t n ) − L ( t n − ) = − M X j =1 v nh ( ξ j ) · Z t n t n − [ τ h ] j d t = − Z t n t n − Z I κ nh · v nh dΓ h ( t )d t, ψ h = v mh into (10c). Bysetting φ h = p nh in (10a) and ϕ h = q nh in (10b), we have Z t n t n − Z I κ nh · v nh dΓ h ( t )d t = Z t n t n − Z I | ∂ Γ ( κ nh · ν h ) | dΓ h ( t )d t. Finally, by combining the two equalities above, we obtain L ( t n ) − L ( t n − ) = − Z t n t n − Z I | ∂ Γ ( κ nh · ν h ) | dΓ h ( t )d t ≤ . This proves that the perimeter of the computed curve decreases in time. (cid:3)
The nonlinear system (14) can be solved by Newton’s iteration. We de-note by ( X nh,l − , p nh,l − , q nh,l − ) the numerical solution obtained in the ( l − th iteration, and define X h,l − ( ξ, t ) = t n − tτ X n − h + t − t n − τ X nh,l − (15) τ h,l − = ∂ ξ X h,l − | ∂ ξ X h,l − | for t ∈ [ t n − , t n ] . (16)In the l th iteration, we calculate ν nh,l − = Z t n t n − ( ∂ ξ X h,l − ) ⊥ d t (17) τ nh,l − = Z t n t n − ∂ ξ X h,l − d t (18) ρ nh,l − = Z t n t n − | ∂ ξ X h,l − | − d t (19) θ nj,l − = Z t n t n − [ τ h,l − ] j d t (20) γ nh,l − = Z t n t n − | ∂ ξ X h,l − | − ∂ ξ X h,l − d t (21) a nh,l − = Z t n t n − | ∂ ξ X h,l − | − t − t n − τ d t (22) A nh,l − = Z t n t n − | ∂ ξ X h,l − | − ( ∂ ξ X h,l − ) T ( ∂ ξ X h,l − ) d t (23)9nd solve the following linear system to obtain the Newton direction ( X δ , p δ , q δ ) ∈ V h × V h × V h , Z X δ τ · (cid:0) φ h ν nh,l − (cid:1) d ξ + Z X nh,l − − X n − h τ · (cid:16) τ φ h ( ∂ ξ X δ ) ⊥ (cid:17) d ξ (24a) − Z ∂ ξ p δ ∂ ξ φ h ρ nh,l − d ξ + Z ∂ ξ p nh,l − ∂ ξ φ h (cid:0) γ nh,l − · ∂ ξ X δ (cid:1) d ξ = − Z X nh,l − − X n − h τ · (cid:0) φ h ν nh,l − (cid:1) d ξ + Z ∂ ξ p nh,l − ∂ ξ φ h ρ nh,l − d ξ ∀ φ h ∈ V h , Z X δ τ · (cid:0) ϕ h τ nh,l − (cid:1) d ξ + Z X nh,l − − X n − h τ · (cid:16) τ ϕ h ∂ ξ X δ (cid:17) d ξ (24b)= − Z X nh,l − − X n − h τ · (cid:0) ϕ h τ nh,l − (cid:1) d ξ ∀ ϕ h ∈ V h , Z p δ ν nh,l − · ψ h d ξ + Z τ p nh,l − ( ∂ ξ X δ ) ⊥ · ψ h d ξ (24c)+ Z q δ τ nh,l − · ψ h d ξ + Z τ q nh,l − ∂ ξ X δ · ψ h d ξ − M X j =1 ψ h ( ξ j ) · [ ∂ ξ X δ a nh,l − ] j + M X j =1 (cid:2) A nh,l − : (cid:0) ψ h ( ξ j ) T ∂ ξ X δ (cid:1)(cid:3) j = − Z (cid:0) p nh,l − ν nh,l − + q nh,l − τ nh,l − (cid:1) · ψ h d ξ + M X j =1 ψ h ( ξ j ) · θ nj,l − ∀ ψ h ∈ V h . The above linearized problem (24) is obtained by using the first-order Taylorexpansion at the point ( X nh,l − , p nh,l − , q nh,l − ) of the nonlinear system (14), andby denoting X δ = X n − X nh,l − , p δ = p n − p nh,l − and q δ = q n − q nh,l − .The algorithm for solving (14) by Newton’s iteration is stated as follows: • Step 1. Choose three proper tolerances TOL , TOL ′ , TOL ′′ >
0. Take theinitial guess of Newton’s iteration to be X nh, = X n − h , p nh, = p n − h and q nh, = q n − h . • Step 2. For given ( X n − h,l − , p n − h,l − , q n − h,l − ), calculate the quantities in (17)-(23) and solve the linear system (24) to obtain ( X δ , p δ , q δ ) ∈ V h × V h × V h .Then set X nh,l = X nh,l − + X δ , p nh,l = p nh,l − + p δ and q nh,l = q nh,l − + q δ . Step 3. If k X nh,l − X nh,l − k L ∞ ≤ TOL, k p nh,l − p nh,l − k L ∞ ≤ TOL ′ and k q nh,l − q nh,l − k L ∞ ≤ TOL ′′ , then set X nh = X nh,l , p nh = p nh,l and q nh = q nh,l and stop; otherwise, set l = l + 1 and goto Step 2.
3. Numerical results
In this section, we present several numerical experiments to demonstratethe accuracy of the proposed numerical scheme (14) and the effectiveness of themethod in preserving the two geometric structures.
We consider the surface diffusion flow of a closed curve with initial shapebeing an ellipse x + 4 y = 4 . We approximate the surface diffusion flow by the proposed numerical scheme(14) and present the relative errors in Table 1, with e h,τ ( t ) = k X h,τ − X h/ ,τ/ k L ∞ ( I ) , (25)where X h,τ is the numerical solution obtained with mesh size h := 1 /M andtime stepsize τ . The order of convergence is calculated byorder of convergence = log (cid:18) k X h,τ − X h/ ,τ/ k L ∞ ( I ) k X h/ ,τ/ − X h/ ,τ/ k L ∞ ( I ) (cid:19)(cid:30) log(2) , based on the finest three meshes. The nonlinear system is solved by usingNewton’s iteration with a tolerance error of 10 − .From Table 1 we can see that the proposed numerical scheme has second-order convergence in space when τ = O ( h ). This is the same as the semi-implicit parametric FEM considered in [4] (see [4, Table 1 on p. 389]).Figure 1 shows the number of Newton’s iterations in the evolution process,with three different time stepsizes. One can see that the method needs onlya few iterations to attain the desired accuracy, and the time stepsize does nothave much influence on the numerical of iterations.11 able 1: Convergence rates in the L ∞ norm, with h = 1 / τ = 0 .
04. The initial curveis an ellipse x + 4 y = 4. h = 1 / τ = 1 / h = 1 / τ = 1 / h = 1 / τ = 1 / h = 1 / τ = 1 / t N = 0 . t N = 0 . t N = 2 1.41E-2 4.08E-3 1.03E-3 2.57E-4 2.00 It e r a t i on N u m be r Figure 1: The number of Newton’s iterations in each time step.
For the example considered in the last subsection, we present in Figure 2the evolution of the normalized perimeter L ( t n ) /L (0), the normalized enclosedarea A ( t n ) /A (0), and a mesh distribution function Ψ( t n ), defined by [4]Ψ( t = t n ) = max ≤ j ≤ M (cid:12)(cid:12) x nj − x nj − (cid:12)(cid:12) min ≤ j ≤ M (cid:12)(cid:12) x nj − x nj − (cid:12)(cid:12) . From Figure 2, one can see that the proposed method can preserve theconservation of A ( t ) /A (0) and the decrease of L ( t ) /L (0), while the mesh distri-bution function Ψ( t ) grows from 1 to around 2 .
3. Figure 3 (a) shows that thearea enclosed by the curve is conserved with machine precision by the proposedmethod; Figure 3 (b) shows that the parametric FEM considered in [4] andBarrett et al. [5] does not strictly preserve the area conservation.12 (a) (b)
Figure 2: Numerical solution given by the proposed method with h = 1 /
32 and τ = 10 − .(a) The normalized perimeter and normalized area enclosed by the curve;(b) The mesh distribution function Ψ( t ). -15 (a) proposed method -5 (b) parametric FEM Figure 3: Comparison between the proposed method and the parametric FEM in [4].Both methods use h = 1 /
32 and τ = 10 − . .3. Evolution of a “flower” shape We apply our method to the numerical simulation of surface diffusion flowwhich initially is a “flower” shape, given in the polar coordinate by x = [1 + 0 .
65 sin(7 θ )] cos θ,y = [1 + 0 .
65 sin(7 θ )] sin θ, θ ∈ [0 , π ] . (26)The shapes of the curves are presented in Figure 4 at several different timelevels. The evolution of the normalized area, normalized perimeter, and meshdistribution function is presented in Figure 5. The two figures show that theflower-shape curve evolves gradually to a circle with the same area and shorterperimeter. -1 0 1-101 (a) -1 0 1-101 (b) -1 0 1-101 (c)-1 0 1-101 (d) -1 0 1-101 (e) -1 0 1-101 (f) Figure 4: Several steps in the evolution of an initially “flower” shape toward its equilibriumat different times: (a) t = 0; (b) t = 10 − ; (c) t = 0 . t = 0 . t = 0 .
006 and (f) t = 0 .
01, where M = 210 and τ = 10 − .
4. Conclusion
We have proposed a novel numerical method for simulating the isotropicsurface diffusion flow of closed curves in two dimensions. Based on a new weak14 (a) (b)
Figure 5: The corresponding temporal evolutions shown in Fig. 4 for: (a) the normalized totalfree energy and the normalized area; (b) the mesh distribution function Ψ( t ), respectively. formulation of the surface diffusion equation, we have introduced a piecewiselinear finite element discretization for the weak formulation to obtain a non-linear system of equations, which can be solved by Newton’s iteration. Wehave rigorously proved that the proposed numerical method can simultaneouslypreserve the two geometric structures, i.e., area conservation and perimeter de-crease. Through numerical experiments it is shown that the proposed method issecond-order convergent in the L ∞ -norm, and can retain the area conservationand perimeter decrease with machine precision. For a complex flower-shapecurve, we have obtained satisfactory numerical results by combining the pro-posed method with a mesh redistribution technique. References [1] Eberhard B¨ansch, Pedro Morin, and Ricardo H. Nochetto. Surface diffusionof graphs: Variational formulation, error analysis, and simulation.
SIAMJ. Numer. Anal. , 42(2):773–799, 2004.[2] Eberhard B¨ansch, Pedro Morin, and Ricardo H. Nochetto. A finite ele-ment method for surface diffusion: the parametric case.
J. Comput. Phys. ,203:321–343, 2005. 153] Weizhu Bao, Wei Jiang, David J Srolovitz, and Yan Wang. Stable equi-libria of anisotropic particles on substrates: a generalized Winterbottomconstruction.
SIAM J. Appl. Math. , 77(6):2093–2118, 2017.[4] Weizhu Bao, Wei Jiang, Yan Wang, and Quan Zhao. A parametric finiteelement method for solid-state dewetting problems with anisotropic surfaceenergies.
J. Comput. Phys. , 330:380–400, 2017.[5] John W. Barrett, Harald Garcke, and Robert N¨urnberg. A parametricfinite element method for fourth order geometric evolution equations.
J.Comput. Phys. , 222:441–467, 2007.[6] John W. Barrett, Harald Garcke, and Robert N¨urnberg. Numerical approx-imation of anisotropic geometric evolution equations in the plane.
IMA J.Numer. Anal. , 28(2):292–330, 2008.[7] John W. Barrett, Harald Garcke, and Robert N¨urnberg. On the parametricfinite element approximation of evolving hypersurfaces in R . J. Comput.Phys , 227(9):4281–4307, 2008.[8] John W. Barrett, Harald Garcke, and Robert N¨urnberg. A variational for-mulation of anisotropic geometric evolution equations in higher dimensions.
Numer. Math. , 109(1):1–44, 2008.[9] John W. Barrett, Harald Garcke, and Robert N¨urnberg. The approximationof planar curve evolutions by stable fully implicit finite element schemesthat equidistribute.
Numer. Methods Partial Differ. Equ. , 27(1):1–30, 2011.[10] John W. Barrett, Harald Garcke, and Robert N¨urnberg. Finite elementmethods for fourth order axisymmetric geometric evolution equations.
J.Comp. Phys. , 376:733–766, 2019.[11] John W. Barrett, Harald Garcke, and Robert Nurnberg. Parametric finiteelement approximations of curvature driven interface evolutions.
Handbookof Numerical Analysis , 21:275–423, 2020.1612] John W. Cahn and Jean E. Taylor. Surface motion by surface diffusion.
Acta Metall. Mater. , 42(4):1045–1063, 1994.[13] B. D. Coleman, R. S. Falk, and M. Moakher. Space–time finite elementmethods for surface diffusion with applications to the theory of the stabilityof cylinders.
SIAM J. Sci. Comput. , 17:1434–1448, 1996.[14] Klaus Deckelnick, Gerhard Dziuk, and Charles M. Elliott. Error analysis ofa semidiscrete numerical scheme for diffusion in axially symmetric surfaces.
SIAM J. Numer. Anal. , 41(6):2161–2179, 2003.[15] Klaus Deckelnick, Gerhard Dziuk, and Charles M. Elliott. Fully discretefinite element approximation for anisotropic surface diffusion of graphs.
SIAM J. Numer. Anal. , 43(3):1112–1138, 2005.[16] K. Deckelnick, G. Dziuk, and C. M. Elliott. Computation of geometricpartial differential equations and mean curvature flow.
Acta Numerica ,14:139–232, 2005.[17] Gerhard Dziuk. An algorithm for evolutionary surfaces.
Numer. Math. ,58(1):603–611, 1990.[18] Gerhard Dziuk. Discrete anisotropic curve shortening flow.
SIAM J. Nu-mer. Anal. , 36(6):1808–1830, 1999.[19] G. Dziuk and C. M. Elliott. Finite elements on evolving surfaces,
IMA J.Numer. Anal. , 27:262–292, 2007.[20] Morton E. Gurtin and Michel Jabbour. Interface evolution in three dimen-sions with curvature-dependent energy and surface diffusion: Interface-controlled evolution, phase transitions, epitaxial growth of elastic films.
Arch. Ration. Mech. Anal. , 163(3):171–208, 2002.[21] Wei Jiang, Weizhu Bao, Carl V. Thompson, and David J. Srolovitz. Phasefield approach for simulating solid-state dewetting problems.
Acta Mater. ,60(15):5578–5592, 2012. 1722] Wei Jiang, Yan Wang, David J Srolovitz, and Weizhu Bao. Solid-statedewetting on curved substrates.
Phys. Rev. Mater. , 2:113401, 2018.[23] Wei Jiang, Yan Wang, Quan Zhao, David J Srolovitz, and Weizhu Bao.Solid-state dewetting and island morphologies in strongly anisotropic ma-terials.
Scripta Materialia , 115:123–127, 2016.[24] Wei Jiang and Quan Zhao. Sharp-interface approach for simulating solid-state dewetting in two dimensions: A Cahn-Hoffman ξ -vector formulation. Physica D: Nonlinear Phenomena , 390:69–83, 2019.[25] Wei Jiang, Quan Zhao, and Weizhu Bao. Sharp-interface model for sim-ulating solid-state dewetting in three dimensions.
SIAM J. Appl. Math. ,80(4):1654–1677, 2020.[26] Bo Li, John Lowengrub, Andreas R¨atz, and Axel Voigt. Geometric evo-lution laws for thin crystalline films: modeling and numerics.
Commun.Comput. Phys. , 6(3):433–482, 2009.[27] William W. Mullins. Theory of thermal grooving.
J. Appl. Phys. ,28(3):333–339, 1957.[28] Jean E. Taylor and John W. Cahn. Linking anisotropic sharp and diffusesurface motion laws via gradient flows.
J. Stat. Phys. , 77:183–197, 1994.[29] Carl V. Thompson. Solid-state dewetting of thin films.
Annu. Rev. Mater.Res. , 42:399–434, 2012.[30] Solmaz Torabi, John Lowengrub, Axel Voigt, and Steven Wise. A newphase-field model for strongly anisotropic systems.
Proc. R. Soc. A ,465:1337–1359, 2010.[31] Yan Wang, Wei Jiang, Weizhu Bao, and David J. Srolovitz. Sharp interfacemodel for solid-state dewetting problems with weakly anisotropic surfaceenergies.
Phys. Rev. B , 91:045303, Jan 2015.1832] Quan Zhao, Wei Jiang, and Weizhu Bao. A parametric finite elementmethod for solid-state dewetting problems in three dimensions.