Ruihan Guo

Francis Filbet

A p-ADAPTIVE LOCAL DISCONTINUOUS GALERKIN LEVEL SET METHOD FOR WILLMORE FLOW

Keywords: Willmore flow, local discontinuous Galerkin method, p-adaptive, semi-implicit scheme, level set method

The level set method is often used to capture interface behavior in two or three dimensions. In this paper, we present a combination of local discontinuous Galerkin (LDG) method and level set method for simulating Willmore flow. The LDG scheme is energy stable and mass conservative, which are good properties comparing with other numerical methods.

In addition, to enhance the efficiency of the proposed LDG scheme and level set method, we employ a p-adaptive local discontinuous Galerkin technique, which applies high order polynomial approximations around the zero level set and low order ones away from the zero level set.

A major advantage of the level set method is that the topological changes are well defined and easily performed. In particular, given the stiffness of Willmore flow, a high order semi-implicit Runge-Kutta method is employed for time discretization, which allows larger time step. The equations at the implicit time level are linear, we demonstrate an efficient and practical multigrid solver to solve the equations. Numerical examples are given to illustrate the combination of the LDG scheme and level set method provides an efficient and practical approach when simulating the Willmore flow.

The level set method was first introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] in 1987, and it has had large impact on computational methods for interface motion. In the past years, the level set method attracted a lot of attention and was used for a wide variety of problems, including fluid dynamics, material sciences, computer vision and imaging processing. For an extensive exposition of the level set method, we refer readers to the review paper [START_REF] Osher | Level set methods: An overview and some recent results[END_REF] and the book [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF].

In the level set method, the evolving family of interface Γ t , t ≥ 0 is represented implicitly as the zero level set of a continuous function which we will denote as φ : Ω × [0, T ] → R, where Ω ⊂ R 2 is a simply connected domain containing the whole family of evolving interfaces Γ t , t ∈ [0, T ]. Therefore, the interface Γ t can be written as

Γ t = x ∈ R 2 , φ(x, t) = 0 .
There are certain advantages associated with this implicit Eulerian type of representation.

• The level set formulation is unchanged in higher dimensions, that is, it is easy to extend to 3D. • Topological changes such as merging and breaking in the evolving interface Γ t are handled naturally.

The level set formulation for Willmore flow will typically result in a fourth order nonlinear partial differential equation (PDE), e.g. see [START_REF] Droske | A level set formulation for Willmore flow[END_REF] for details. Beneš et al. developed a finite volume method for spatial discretization in [START_REF] Beneš | Comparison study for level set and direct Lagrangian methods for computing Willmore flow of closed planar curves[END_REF]. When the evolved surface is a graph, the Willmore flow converts to Willmore flow of graphs, and numerical methods for Willmore flow of graphs attracted a lot of attention in the last decades. A semi-implicit numerical scheme for the Willmore flow of graphs with continuous finite element discretization and the convergence analysis has been provided by Deckelnick and Dziuk in [START_REF] Deckelnick | Error analysis of a finite element method for the Willmore flow of graphs[END_REF]. Finite difference discretization for the Willmore flow of graphs has been discussed in [START_REF] Oberhuber | Numerical solution for the Willmore flow of graphs[END_REF]. Xu and Shu [START_REF] Xu | Local discontinuous Galerkin method for surface diffusion and Willmore flow of graphs[END_REF] developed a local discontinuous Galerkin method for Willmore flow of graphs, and time discretization was by the forward Euler method with a suitably small time step for stability.

The discontinuous Galerkin (DG) method is a class of finite element methods using completely discontinuous basis functions, which are usually chosen as piecewise polynomials. It was first introduced in 1973 by Reed and Hill [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF] for solving the linear time dependent neutron transport equation. Later, the DG method was developed for solving the nonlinear hyperbolic conservation laws by Cockburn et al. in a series of papers [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one dimensional systems[END_REF][START_REF] Cockburn | The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case[END_REF][START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems[END_REF].

These discontinuous Galerkin methods were generalized to solve the PDEs containing higher than first order derivatives, i.e. the local discontinuous Galerkin (LDG) method. The first LDG method was designed to solve a convection diffusion equation (with second derivatives) by Cockburn and Shu [START_REF] Cockburn | The local discontinuous Galerkin method for time-dependent convectiondiffusion systems[END_REF]. The idea of LDG methods is to suitably rewrite a higher order PDE into a first order system, then apply the DG method to the system. A key ingredient for the success of such methods is the correct design of interface numerical fluxes, which should be designed to guarantee stability. For a detailed description about the LDG methods for high order time dependent PDEs, we refer the readers to [START_REF] Xu | Local discontinuous Galerkin methods for high-order time-dependent partial differential equations[END_REF]. The DG and LDG methods have several attractive properties, for example:

(1) The order of accuracy can be locally determined in each cell, thus allowing for efficient p adaptivity. (2) These methods can be used on arbitrary triangulations, even those with hanging nodes, thus allowing for efficient h adaptivity. (3) These methods have excellent parallel efficiency.

In this paper, we present the level set method for computing evolution of the Willmore flow of embedded surfaces, and one of the major advantages of level set method is their ability to easily handle topological changes. Then we develop a high order local discontinuous Galerkin method for solving level set formulation of the Willmore flow. The basis functions of LDG method can be completely discontinuous, which leads that the order of accuracy can be locally determined in each cell, i.e. p-adaptivity. To improve the efficiency of the proposed methods, a p-adaptive technique is applied here to capture the evolution of the interface.

Another main difficulty when solving the Willmore flow numerically is that it is stiff, and explicit time discretization methods suffer from extremely small time step restriction for stability, of the form ∆t ≤ C∆x 4 , which is not efficient, especially for long time simulation. It would therefore be desirable to develop implicit or semi-implicit time marching techniques to alleviate this problem. In [START_REF] Xia | Efficient time discretization for local discontinuous Galerkin methods[END_REF], three efficient time discretization techniques were explored, these are the semi-implicit spectral deferred correction (SDC) method, the additive Runge-Kutta (ARK) method and the exponential time differencing (ETD) method. These methods are mainly based on a splitting of stiff and non stiff differential operators, which is particularly efficient for semilinear PDEs, but for quasi-linear equations, like for Willmore flows, fully nonlinear solver would be required. Thus we will employ here a high order semi-implicit Runge-Kutta method [START_REF] Boscarino | High order semi-implicit schemes for time-dependent partial differential equations[END_REF][START_REF] Guo | Efficient high order semi-implicit time discretization and local discontinuous Galerkin methods for highly nonlinear PDEs[END_REF].

The organization of the paper is as follows. In Section 2, we review the level set formulation for Willmore flow and get a highly nonlinear fourth order PDE. In Section 3, we present the formulation of the local discontinuous Galerkin method and the p-adaptive technique for solving the Willmore flow. Section 4 contains a simple description of the high order semiimplicit Runge-Kutta method. Numerical experiments are presented in Section 5, testing the performance of p-adaptive local discontinuous Galerkin level set method for Willmore flow. Finally, we give concluding remarks in Section 6.

Level set formulation for Willmore flow

In the level set method, the evolving family of interface Γ t , t ≥ 0 is represented implicitly as the zero level set of a continuous function φ. Therefore, the interface Γ t can be written as

Γ t = x ∈ R 2 , φ(x, t) = 0 .
For the Willmore flow, the level set function

∂φ ∂t + V |∇φ| = 0,
with V being the speed of propagation of the level set function φ, turns into the following equation

(2.1) ∂φ ∂t + ∆ Γt H(t) + H(t) S(t) 2 L 2 - 1 2 H(t) 2 |∇φ| = 0,
where S(t) denotes the shape operator on Γ t , H(t) is the mean curvature on Γ t and is given as

H = ∇ • ∇φ |∇φ| .
In order to obtain the Willmore flow, it is necessary to reformulate the equation in terms of φ(t), H(t) and their derivatives.

For the sake of simplicity, we assume that |∇φ| = 0 on Γ t . Then the unit normal vector to a level curve of φ can be defined as

(2.2) n = ∇φ |∇φ| .
Let us denote the following auxiliary functions Q = |∇φ| and ω = Q H. Then, by the derivation in [START_REF] Droske | A level set formulation for Willmore flow[END_REF], a level set formulation of Willmore flow is given as follows:

(2.3)

           ∂φ ∂t + Q ∇ • E ∇ω - 1 2 ω 2 Q 3 ∇φ = 0, ω = Q ∇ • ∇φ |∇φ| .
where the matrix E is given by

E = 1 Q I - ∇φ Q ⊗ ∇φ Q ,
which is a projection into a tangential space of the curve representing the zero level set of φ.

The Willmore flow (2.3) is subjected to the initial condition

φ(x, 0) = φ 0 (x), x ∈ Ω,
where the initial function φ 0 (x) is a signed distance function, i.e. φ 0 (x) = dist(x, Γ 0 ), for any x ∈ Ω. The initial condition satisfies the following property

(2.4) |∇φ 0 (x)| = 1.
One of the major advantages of level set methods is their ability to easily handle topological changes. However, we have found this not to be the case for Willmore flow because of the lack of a maximum principle. Indeed, two surfaces both undergoing an evolution by Willmore flow may intersect in finite time. Therefore, the level set formulation for Willmore flow will lead to singularities in general, i.e. ∇φ = 0 in finite time. To handle these difficulties, it is necessary to introduce a suitable regularization technique. For a fixed regularization parameter 0 < ε 1, we define (2.5)

Q ε = ε + |∇φ| 2
and replace all occurrences of Q in (2.3) by Q ε .

Numerical method

With the level set formulation for Willmore flow, we need to discrete equation (2.3) to capture the movement of the interface (zero level set). It would therefore be desirable to develop high order numerical scheme in both space and time to obtain high resolution. In this section, we will consider the local discontinuous Galerkin (LDG) method for the Willmore flow (2.3).

3.1. Notation. We consider a subdivision T h of Ω with shape-regular elements K. Let E h denote the union of the boundary faces of elements

K ∈ T h , i.e. E h = ∪ K∈T h ∂K, and E 0 = E h \ ∂Ω.
In order to describe the flux functions, we need to introduce some notations. Let e be an interior face shared by the "left" and "right" elements K L and K R and define the normal vectors ν L and ν R on e pointing exterior to K L and K R , respectively. For our purpose "left" and "right" can be uniquely defined for each face according to any fixed rule. For example, we choose ν 0 as a constant vector. The left element K L to the face e requires that ν L • ν 0 < 0, and the right one K R requires ν L • ν 0 ≥ 0. If ψ is a function on K L and K R , but possibly discontinuous across e, let ψ L denote (ψ| K L )| e and ψ R denote (ψ| K R )| e , the left and right trace, respectively.

Let P k (K) be the space of polynomials of degree at most k ≥ 0 on K ∈ T h . The finite element spaces are denoted by

V h = ϕ ∈ L 2 (Ω), ϕ| K ∈ P k (K), ∀ K ∈ T h and Σ d h = Φ = (φ 1 , • • • , φ d ) T ∈ (L 2 (Ω)) d , φ l | K ∈ P k (K), ∀l = 1 • • • d, ∀ K ∈ T h .
Note that functions in V h , Σ d h are allowed to be completely discontinuous across element interfaces.

Local discontinuous Galerkin method.

To develop the local discontinuous Galerkin scheme, we first rewrite the Willmore flow (2.3) as a first order system:

(3.1)            1 Q ε ∂φ ∂t + ∇ • (s -v) = 0, s = E p and v = 1 2 ω 2 Q 3 ε r,
with p = ∇ω, where ω is given by the algebraic equation ω = Q ε H, and H is solution to H = ∇ • q, where q satisfies Q ε q = r and finally r = ∇φ. In equation (3.1), the matrix E is given by

(3.2)          E = 1 Q ε I - r ⊗ r Q 2 ε , Q ε = ε + |r| 2 .
Next, we obtain the weak formulation of the exact solution. Applying the LDG method to system (3.1), we get the following numerical scheme for the unknowns (φ, ω, H) ∈ V 3 h and (s, v, p, q, r) ∈ Σ 5 h , such that, for all test functions (ϕ, ξ, ϑ) ∈ V 3 h and (θ, ψ, η, ρ, ζ) ∈ Σ 5 h , we have

(3.3) K 1 Q ε ∂φ ∂t ϕ dK = K (s -v) • ∇ϕ dK - ∂K ( s • ν -v • ν) ϕ ds,
where the previous algebraic relations between (s, v) and (p, r, ω) are now written as

           K s • θ dK = K E p • θ dK, K v • ψ dK = K 1 2 ω 2 Q 3 ε r • ψ dK.
Moreover, p is solution to

K p • η dK = - K ω ∇ • η dK + ∂K ων • η ds,
where ω is now given with respect to

H K ωξ dK = K Q ε H ξ dK,
and H is solution to

K Hϑ dK = - K q • ∇ϑ dK + ∂K q • νϑ ds,
where (q, r) are given with respect to the main unknown φ, that is,

           K q • ρ dK = K r Q ε • ρ dK, K r • ζ dK = - K φ∇ • ζ dK + ∂K φν • ζ ds,
where E and Q ε are computed by (3.2). Here s, v, ω, q and φ are the numerical fluxes. To complete the definition of the LDG method, we need to define these numerical fluxes, which are functions defined on the edges and should be designed based on different guiding principles for different PDEs to ensure stability.

Similar to the development in [START_REF] Xu | Local discontinuous Galerkin method for surface diffusion and Willmore flow of graphs[END_REF], it turns out that we can take the simple choices such that (3.4) s| e = s L , v| e = v L , q| e = q R , ω| e = ω L , φ| e = φ R .

We remark that the choice for the numerical fluxes (3.4) is not unique. In fact, the crucial part is taking s, v and φ from opposite sides, and q and ω from opposite sides.

Remark 3.1. The solution to the LDG scheme is mass conservative, by choosing the test function ϕ = 1 in equation (3.3).

Remark 3.2. Similar to the proof in [START_REF] Xu | Local discontinuous Galerkin method for surface diffusion and Willmore flow of graphs[END_REF], the solution to the LDG scheme (3.3) with numerical fluxes (3.4) satisfies the following energy stability 1 2

d dt Ω H 2 Q ε dx + Ω (φ t ) 2 Q ε dx = 0.

A p-adaptive technique.

In the level set method, we aim to capture the interface behavior. As a result, we are mainly interested in the position near the zero level set, which motivates us to use adaptive technique. The basis functions of the LDG method we discuss in this paper can be discontinuous, so it has certain adaptive flexibility, such as, the order of accuracy can be locally determined in each cell.

Our aim is to apply a higher order polynomial approximation around the zero level set to achieve high order accuracy, whereas in the regions far away from the zero level set, we apply lower order polynomial approximations, which is similar as what is described by Yan and Osher in [START_REF] Yan | Discontinuous Galerkin level set method for interface capturing[END_REF]. We explain the p-adaptive technique in detail with Figure 3.1. In Figure 3.1, we apply P 2 approximation near the zero level set, which are the domain filled with dark cells. While for the domain two or three cells away from the zero level set, it is not necessary to use high order approximation, thus P 1 approximation is enough.

As we know, in the implementation of the p-adaptive technique, it is essential to employ a flag to mark "near the zero level set". By the definition of the level set method, it is reasonable to use the following flag:

P 2 approximation, if |φ| ≤ 1.5h, P 1 approximation, if |φ| > 1.5h, (3.5)
where |φ| is the absolute value of φ. With the proposed p-adaptive technique, we can improve the efficiency of the local discontinuous Galerkin method for level set equation.

Time discretization

As we know, the Willmore flow (2.3) is a highly nonlinear fourth order PDE, explicit time integration methods lead to an extremely small time step restriction for stability, but not for accuracy. It would therefore be desirable to develop a semi-implicit time discretization technique to alleviate this problem. In this section, we employ a high order semi-implicit Runge-Kutta method. Based on the work in [START_REF] Boscarino | High order semi-implicit schemes for time-dependent partial differential equations[END_REF], the main idea of the semi-implicit method is to apply two different Runge-Kutta methods. Consider the ODE system (4.1)

   du dt = H(t, u(t), u(t)), t ∈ [0, T ] u(0) = u 0 ,
where the function H : R × R m × R m → R m is sufficiently differentiable and the dependence on the second argument of H is non-stiff, while the dependence on the third argument is stiff. The semi-implicit Runge-Kutta method is based on the partitioned Runge-Kutta method, so we will first give a simple description for the partitioned Runge-Kutta methods. We consider autonomous differential equations in the partitioned form, (4.2)

     dy dt (t) = F(y(t), z(t)), dz dt (t) = G(y(t), z(t)),
where For practical reasons, in order to simplify the computations, we consider that  is a strictly lower triangular matrix and A is a lower triangular matrix for the implicit part. In addition, the coefficients satisfy:

y(t) ∈ R m , z(t) ∈ R n , F : R m × R n → R m , G : R m × R n → R n and F, G ∈ C 1 (R m × R n ). y(t 0 ) = y 0 , z(t 0 ) =
(4.3) ĉi = i-1 j=1
âi,j , and c i = i j=1 a i,j , for 1 i s.

By applying the partitioned Runge-Kutta time marching method, the solution of the autonomous system (4.2) advanced from time t n to t n+1 = t n + ∆t is given by (4.4)

             k i = F y n + ∆t i-1 j=1 âi,j k j , z n + ∆t i j=1 a i,j l j , 1 i s, l i = G y n + ∆t i-1 j=1 âi,j k j , z n + ∆t i j=1 a i,j l j , 1 i s,
and we can calculate y n+1 and z n+1 as follows (4.5)

           y n+1 = y n + ∆t s i=1 bi k i , z n+1 = z n + ∆t s i=1 b i l i .
To derive a semi-implicit Runge-Kutta scheme, we first rewrite the non autonomous differential equation (4.1) as an autonomous system where we double the number of variable, that is,

(4.6)                d dt t u(t) = 1 H(t, u(t), v(t)) , dv(t) dt = H(t, u(t), v(t)), u(0) = v(0) = u 0
By a uniqueness argument the solution to (4.6) corresponds to the one of (4.1), but now this system is written as an autonomous partitioned system (4.2), with y(t) = (t, u(t)), F = (1, H) and z(t) = u(t), G = H and y(t 0 ) = (t 0 , u 0 ), z(t 0 ) = u 0 are the initial conditions. Applying the partitioned Runge-Kutta scheme (4.4)-(4.5) to system (4.6), we can get a high order semi-implicit Runge-Kutta method for (4.

k i = H t n + ĉi ∆t, u n + ∆t i-1 j=1 âi,j k j , u n + ∆t i j=1 a i,j k j , 1 i s.
It is worth to mention here that k i is defined implicitly since a i,i = 0. Therefore, starting from u n , we give the algorithm to calculate u n+1 in the following.

(1) Set for i = 1, . . . , s,

U i = u n + ∆t i-1 j=1 âi,j k j , V i = u n + ∆t i j=1 a i,j k j . (4.7) (2) For i = 1, . . . , s, compute (4.8) k i = H(t n + ĉi ∆t, U i , V i ).
(3) Update the numerical solution u n+1 as (4.9)

u n+1 = u n + ∆t s i=1 b i k i .
After the LDG spatial discretization for the Willmore flow (2.3), we can apply the semiimplicit scheme (4.7)-(4.9) by writing the system of ODE in the form of (4.1) with the component u(t) treated explicitly, the component v(t) treated implicitly and

   H(t, u, v) = -Q ε (u) ∇ • E(u)∇ω(u, v) -1 2 ω(u) 2 Qε(u) 3 ∇v , ω(u, v) = Q ε (u) ∇ • ∇v Qε(u) , where    Q ε (u) = ε + |∇u| 2 , ω(u) = Q ε (u) ∇ • ∇u Qε(u) , E(u) = 1 Qε(u) I -∇u Qε(u) ⊗ ∇u Qε(u) .
By algorithm (4.7)-(4.9), the second variable of equation (4.1) is treated explicitly and the third one is treated implicitly. Obviously, it is necessary to solve system of linear algebraic equations (4.8) at each time step. The overall performance highly depends on the performance of the solver. Traditional iterative solution methods such as Gauss-Seidel method suffers from slow convergence rates, especially for large systems. To enhance the efficiency of the high order semi-implicit time marching method, we will apply the multigrid solver to solve algebraic equations (4.8) in this paper. And for a detailed description of the multigrid solver, we refer the readers to Trottenberg et al. [START_REF] Trottenberg | Multigrid[END_REF] and Brandt [START_REF] Brandt | Multigrid techniques: 1984 guide with applications to fluid dynamics[END_REF].

Numerical results

In this section, we present some numerical results for the Willmore flow. In addition, we use p-adaptive LDG spatial discretization (P 2 approximation near the zero level set, P 1 approximation away from the zero level set) coupled with a high order semi-implicit Runge-Kutta method, which is high order accurate in both time and space, allowing larger time step. Each time step we solve the linear algebraic equations by multigrid solver. We first present various representative examples of shape relaxation, for example, the motions of an ellipse, a square, an asteroid and a singular curve. Next, we will show the topological changes, such as merging and breaking. 5.1. Accuracy test. In this example, we consider the accuracy test for one-dimensional Willmore flow. We consider the exact solution (5.1) φ(x, t) = 0.05 sin(x) cos(t)

for equation (2.3) with a source term f (x, t), which is a given function so that (5.1) is the exact solution.

To show the effect of different regularization parameters ε on the corresponding errors and orders of accuracy, we report the L 2 and L ∞ errors and the numerical orders of accuracy at time T = 0.5 with uniform meshes in Table 5.1. We can see that the method with P k elements gives (k + 1)-th order of accuracy with different scales of the regularization parameters ε.

ε = 1 ε = 5h ε = 5h 2
Mesh L 2 error order L ∞ error order L 2 error order L ∞ error order L 2 error order L ∞ error order The motion of an ellipse. We consider the motion of an ellipse as shown in Figure 5.1. For large time, the level set asymptotically approaches a circle. Here, we approximate the Willmore flow with ε = h on a 64 × 64 grid for the computational domain Ω = [0, 4] × [0, 4]. Here, we adopt a p-adaptive LDG method for spatial discretization and a second order semi-implicit Runge-Kutta method with a time step of ∆t = 0.001h, which is very large compare to that of an explicit time marching method (∆t = Ch 4 ). The motion of a square. We now consider an initial data, where the initial level set is a square, that is, it has sharp corners. We report the numerical results in Figure 5.2, which shows the evolution of a square into a circle. We approximate the Willmore flow with ε = h on a 64 × 64 grid for the computational domain Ω = [0, 4] × [0, 4]. Here, we also adopt a second order semi-implicit Runge-Kutta method with a time step of ∆t = 0.001h. The motion of an asteroid. In Figure 5.3, we present the motion of another non-convex curve with very sharp corners. The numerical methods and parameters are chosen to be the same with those in the previous example. Also the initial asteroid is evolved to be a circle.

P 1 16
The motion of a singular curve. We consider an initial data, where the initial level set is a singular curve. We report the numerical results in Figure 5.4, which shows the evolution of a singular curve into a circle. We approximate the Willmore flow with ε = 10h and ε = 10h both undergoing an evolution by Willmore flow my intersect in finite time. To handle these difficulties, it is important to choose the regularization parameter ε to avoid a blow up of the gradient of φ in finite time. In addition, it is important to choose small time steps to temporally resolve the rapid transients that occur at a topological change. Figure 5.5 presents the merging of two square-like shapes under the level set evolution of Willmore flow. We can see the two square-like shapes merge and eventually evolve to a single circle. Here, we choose ε = 5h, and the time step size is ∆t = 10 -6 . Since for such small time step, we use here a first order semi-implicit time marching method. Certainly, we can use explicit Runge-Kutta method here, but because of the severe time step restriction of explicit methods, the time step must be very small, which takes about ∆t ≈ 10 -9 .

In Figure 5.6, we show another example of topological change. The initial configuration consists of two almost touching surfaces-the inner surface being a circle and the outer surface being an ellipse. We can see a phenomenon of pinching off in Figure 5.6. Here, we also choose ε = 5h, and the time step size ∆t = 10 -6 .

Concluding remarks

In this work, we have developed a local discontinuous Galerkin method for a level set formulation of Willmore flow. The scheme is mass conservative and energy stable. In the level set method, we are mainly interested in capturing the movement of the interface. Therefore, it is reasonable and natural to combine our local discontinuous Galerkin scheme with p-adaptive technique to improve the efficiency of the method. To relax the severe time step restriction of explicit time integration methods for stability and achieve high order temporal accuracy, a semi-implicit Runge-Kutta method is adopted. The equations at the implicit time level are linear and we employ an efficient multigrid solver to solve them. Numerical experiments demonstrate that the proposed methods are efficient for Willmore flow, and the level set method can easily handle topological changes. (c) t = 0.00022 (f) t = 0.0003

Figure 3 . 1 .

 31 Figure 3.1. Explanation of p-adaptive technique.

  z 0 are the initial conditions. Then we can express the partitioned Runge-Kutta methods by applying two different Runge-Kutta methods as the following

  1) : the first component of the first equation (4.6) only gives ĉi = i j=1 âi,j , whereas the second component of the first equation and the second equation of (4.6) are identical, which gives the following semi-implicit scheme

1 Figure 5 . 1 .

 151 Figure 5.1. Towards a circle. Evolution of Willmore of the regularized flow with ε = h when the initial level set corresponds to an ellipse.

Figure 5 . 2 .

 52 Figure 5.2. Towards a circle. Evolution of Willmore of the regularized flow with ε = h when the initial level set corresponds to a square.

Figure 5 . 3 .

 53 Figure 5.3. Towards a circle. Evolution of Willmore of the regularized flow with ε = h when the initial level set corresponds to an asteroid.

1 Figure 5 . 4 .

 154 Figure 5.4. Towards a circle. Evolution of Willmore of the regularized flow with ε = 10h and ε = 10h 2 when the initial level set corresponds to a singular curve.

Figure 5 . 5 .

 55 Figure 5.5. Topology changes. Two square-like shapes merge under the level set evolution of Willmore flow.

Figure 5 . 6 .

 56 Figure 5.6. Topology changes. Two shapes pinch off under the level set evolution of Willmore flow.

Table 5 .

 5 1. Accuracy test. L 2 and L ∞ error norms on the solution at time T = 0.5 with different regularization parameters ε.5.2. Towards a circle.We first perform a series of numerical simulations, where the initial level sets describe two convex curves, a non-convex curve and a singular curve. Hence it is expected that the level set approaches a circle for large time, which compares very well with numerical calculations performed by Beneš et al.[START_REF] Beneš | Comparison study for level set and direct Lagrangian methods for computing Willmore flow of closed planar curves[END_REF].