AREA PRESERVING GEODESIC CURVATURE DRIVEN FLOW OF CLOSED CURVES ON A SURFACE

. We investigate a non-local geometric ﬂow preserving surface area enclosed by a curve on a given surface evolved in the normal direction by the geodesic curvature and the external force. We show how such a ﬂow of surface curves can be projected into a ﬂow of planar curves with the non-local normal velocity. We prove that the surface area preserving ﬂow decreases the length of the evolved surface curves. Local existence and continuation of classical smooth solutions to the governing system of partial diﬀerential equations is analysed as well. Furthermore, we propose a numerical method of ﬂowing ﬁnite volume for spatial discretization in combination with the Runge–Kutta method for solving the resulting system. Several computational examples demonstrate variety of evolution of surface curves and the order of convergence.

1. Introduction and model description. This article introduces the evolution of the family {G t } t≥0 of closed, non-selfintersecting curves evolving on a given two dimensional surface M ⊂ R 3 . It is assumed that the surface M is represented by the graph of a smooth function ϕ : R 2 → R defined in a domain Ω ⊂ R 2 , i.e., The family of curves {G t } t≥0 evolves according to the following geometric law: where G ini ⊂ M is the initial curve, K G is the geodesic curvature and F is the prescribed external scalar force, L(G t ) = G dS is the length of a surface curve G t .
A particular choice of the forcing term F in (1) corresponds to the area-preserving curvature flow for planar curves Γ t , t ≥ 0, in which the normal velocity has the form: It is known that the flow (3) preserves the area A(Γ t ) enclosed by the curve Γ t , i.e., A(Γ t ) = A(Γ ini ) for all t ≥ 0 (c.f. Gage [14]). In what follows, we show that the analogous conservation property holds for the flow (1), i.e., the area A(G t ) enclosed by the curve G t on the surface M is preserved for all t ≥ 0. Geometric evolution equations of the form (3) [14,11,20,13,17,9]). Recently, in [16] Kolář, Beneš andŠevčovič investigated the geometric law (3) for motion of open planar curves with fixed ends. In the context of the modified Allen-Cahn equation (c.f. [8,1]) approximating the curvature driven flow (3) (see [3]), the constrained curvature driven flow was also studied in [25,7,15,6] by Rubinstein  Motivation for studying the area-preserving constrained motion (3) driven by the mean curvature has its origin in physics of phase transitions, e.g. in recrystallization, where a previously melted, fixed volume of the liquid phase solidifies again (see [19]). It is natural to generalize the area preserving flow of planar curves to the case when the curves are evolved on a given two dimensional surface.
In this paper, problem (1) for closed curves is mathematically treated by the parametric (or also called direct or Lagrangian) method, which has been applied for planar curve motion by various authors. We refer the reader to papers [10,4,5] by Deckelnick and Beneš et al. where the basics of this approach and its applications are elaborated. In [24] Pauš et al. compared the direct approach and other interface-capturing methods, such as the level-set method or the phase-field method. In the paper [24] by Pauš et al., the approximate algorithmic approach for handling topological changes (like self-intersecting) was proposed and analyzed. Concerning some of the drawbacks of the direct method, the problem of tangential redistribution was analyzed in [27] byŠevčovič and Yazaki. The application of the direct approach to the flow (1) results into a system of degenerate parabolic equations for the parametrization of the curve G t . The system of governing equations is solved numerically to provide an information about the behavior of a solution of (1). The numerical approximation scheme is based on the flowing finite volume method which was proposed by Mikula andŠevčovič in [21] for curvature driven flows of planar curves.
The paper is organized as follows. Next section recalls the direct Lagrangian approach for solving curvature driven flows of curves. In Section 3 we derive a system of nonlocal partial differential equations for parametrization of evolving curves. Next in Section 4 we briefly discuss the role of the tangential velocity and we present a class of the curvature adjusted tangential velocities. Local existence, uniqueness and continuation of Hölder smooth classical solutions is shown in Section 5. Section 6 is devoted to the area preserving flow of surface curves. We show that the enclosed area of a family of evolving curves driven in normal direction by the velocity (1) is preserved and their length is shortened. In Section 7 we use the numerical method of flowing finite volumes combined with the Runge-Kutta method for solving governing PDEs. Finally, in Section 8 we present several computational examples of area preserving flow of curves evolving on a given surface.
2. Projection method. In the direct parametric (or Lagrangian) approach, a planar time-dependent curve Γ t ⊂ R 2 , t ≥ 0, is described by the position vector: where u is the parameter from a fixed interval. Since we are concerned with closed curves for which X is 1-periodic in the u variable, we will identify the interval [0, 1] with R/Z ≡ S 1 . Then the curve Γ t is given by: where X 1 (u, t) and X 2 (u, t) are the components of the position vector X(u, t). Now let us consider a closed curve G t on a surface M, which is a graph of a function ϕ : Ω ⊂ R 2 → R. Such a curve can be uniquely represented by its vertical projection to the plane Ω ⊂ R 2 , i.e., where Γ t is a planar curve in Ω ⊂ R 2 (see Figure 1). Our approach is to analyze the flow of curves G t ∈ M on a surface driven by (1) by means of the flow of projected curves Γ t ∈ R 2 in the plane.
In what follows, we will derive a system of governing equations for parametrization X(u, t) of Γ t provided that G t is evolved in the normal direction by (1). To this end, we have to find the normal velocity v Γ of the curve Γ t in terms of geometric quantities corresponding to Γ t . Assume the parametrization of a closed curve Γ t is oriented anticlockwise, and the periodic boundary conditions at u = 0 and u = 1 are imposed, i.e., X| u=0 = X| u=1 and ∂ u X| u=0 = ∂ u X| u=1 . Then, the geometric quantities such as the unit tangential vector t Γ , the outer unit normal vector n Γ , and the curvature κ Γ can be expressed in terms X as follows: Notice that the choice of the normal vector in the outer direction is in accordance with the rule det(n Γ , t Γ ) = 1. Let us denote the arc-length variable by s. Then ds = |∂ u X(u, t)| du. The Frénet formulas yield where the curvature κ Γ is given by the inner product: Here a · b stands for the Euclidean inner product in R 2 . Having a surface curve G t ⊂ M, we consider the important geometrical quantities, such as the unit tangent and outer normal vectors T , N belonging to the tangent space T x (M) to G t ⊂ M, and the geodesic curvature K G . Our aim is to express them in terms of the quantities: n Γ , t Γ , κ Γ , and ∇ϕ. If the surface M = {(x, ϕ(x)) ∈ R 3 : x ∈ Ω} is a graph of the function ϕ then the unit normal vector to M is given by: The vector N M together with the unit tangent T and outer normal N to G t define the moving orthonormal frame. We can express the vectors T and N as follows: Finally, the geodesic curvature K G of the curve G on the surface M is given by where κ Γ is the curvature of the projected planar curve Γ t (see [22,Eq. (6)] where the authors considered the inner normal vector instead).
3. Planar motion law. Our aim is to derive a system of equations for the position vector X(u, t) provided that the corresponding family of surface curves G t satisfies the geometric equation (1). We seek the geometric equation for the normal velocity v Γ of the planar curve Γ t in the form such that Γ t is the vertical projection of G t . For description of the time evolution of the position vector X(u, t) ∈ Γ t , we consider the following geometric equation: where β = β(X, n Γ , κ Γ ) and α are the normal and the tangential components of the velocity of the planar curve Γ t . Notice that the presence of a tangential velocity does not change the shape of evolving closed curves. But it has a strong impact on redistribution of points along the curve (see, e.g., [21]).
The normal velocity V G of a surface curve G t is a projection of the speed of the position vector (X(u, t), ϕ(X(u, t))) T on G t into the normal direction N , i.e., (see [22]). It follows from (1) and (7), that the velocity β is given by: Since the arc-length parametrization S of the surface curve G t satisfies dS = 1 + (∇ϕ · t Γ ) 2 ds, we obtain and K G is given by (7) and L(G t ) = Γ 1 + (∇ϕ · t Γ ) 2 ds. Now it follows from equations (6), (7) and (11) that the curve G t evolves on the surface M according to law (1) provided that the vertically projected planar curve Γ t satisfies the geometric equation (8) in the form: where the coefficients a > 0, b and c are smooth functions defined as The curve Γ t then evolves according to the geometric evolution law (8) provided that its parametrization X(u, t) satisfies the following system of degenerate parabolic equations: where the coefficients a, b and c are given by (12), (13) and (14), and the non-locally defined force term F, is expressed as follows: 4. Tangential velocity. We will briefly discuss the choice of a suitable tangential velocity functional α. To construct a stable numerical computational scheme, several nontrivial choices of α have been proposed in the literature. We refer the reader to papers by Mikula andŠevčovič [22,23], Beneš et al. [5,16,24] and references therein. A general framework yielding the so-called curvature adjusted tangential velocity has been proposed byŠevčovič and Yazaki in [27]. This approach takes into account variations in the curvature as well as the necessity of uniform or asymptotically uniform redistribution of points along evolved curves. More precisely, we shall define the so-called curvature adjusted relative local length quantity: and F = (1/L(Γ t )) Γ F (s) ds is the arc-length average of a quantity F over the curve Γ t . Here Φ : R → R is a suitable nonnegative shape function depending on the curvature. In [27] it was shown that lim t→Tmax r(u, t) = 1 (19) uniformly with respect to u ∈ [0, 1] provided that the tangential velocity α is a solution to the following equation: where Φ = Φ(κ Γ ) and k 1 , k 2 > 0 are some constants. To construct a unique solution α, we assume the renormalization condition Φ(κ Γ )α = 0, (c.f. [28,Eq. (10)]).
The following lemma deals with properties of the tangential velocity functional and it is due toŠevčovič and Yazaki [28]. To formulate its statement we need to introduce the scale of Banach spaces Here c 2k+ε (S 1 ) is the so-called little Hölder space (c.f. Angenent [2]). It is the closure of C ∞ smooth 1-periodic functions in the norm of the Hölder space C 2k+ε for some positive 0 < ε < 1. By c 2k+ε * we denoted the space c 2k+ε * be the tangential velocity function given as a unique solution to (20) satisfying the renormalization condition

5.
Local existence and continuation of classical solutions. In this section we prove local existence, uniqueness and continuation of classical Hölder smooth solutions to the system of governing PDEs (21) below. Throughout this section we shall assume that the mapping φ : Ω ≡ R 2 → R is at least C 5 smooth and it has bounded derivatives up to the second order, First, we prove a useful lemma giving us an a-priori estimate of the external force F term (1).
, be a flow of surface curves with the normal velocity V G given by (1). Then there exists a constant C 1 > 0 depending on the area A(G ini ) enclosed by the initial curve G ini only, and such that Proof. First we recall the Gauss-Bonnet formula: which is satisfied by any closed non-selfintersecting curve G on a simple surface M where K is the Gaussian curvature of M (c.f. [26,Chapter 2]). Hence, where K max = max M |K| is the maximum of the modulus of the Gaussian curvature of the surface M which is bounded due to the assumptions made on the function φ. Since then, by using the isoperimetric inequality L 2 ≥ 4πA in the plane R 2 , we obtain the lower bound for the length L(G t ): Hence, there exists a constant C 1 > 0 depending on the initial area A(G ini ) only and such that 1 as claimed.
We can state the following local existence and uniqueness result.
Proposition 1. Assume the parametrization X ini of an initial surface curve G ini is at least C 4+ε smooth. Suppose that the tangential velocity α belongs to the class of curvature adjusted tangential velocities introduced in Section 4 or α = 0. Then there exists T > 0 and the unique family of surface curves G t , t ∈ [0, T ], evolving in the normal direction with the velocity V G given by (1) and such that its parametriza- If the maximal time of existence is finite, T max < +∞, then lim sup t→Tmax max |K Gt (G t )| = +∞.
Proof. Recall that the normal velocity β(X, n Γ , κ Γ ) of the projected planar curve Γ t in the outer normal direction n Γ can be written in the form β = −aκ Γ + b + cF.
The coefficients a, b, c depends on the position vector X and the unit outer normal n Γ and the tangent vector t Γ , where n Γ = t ⊥ Γ (see (12), (13) and (14)). In [28, Theorem 1]Ševčovič and Yazaki proved a rather general result on local existence and uniqueness and continuation of classical Hölder smooth solutions to the non-local flow driven by the normal velocity which is the sum of local and nonlocal parts provided that the nonlocal part is however independent of X and n Γ . This is why we have to slightly modify the proof of [28,Theorem 1] in order to handle the case when the nonlocal part has the form: c(X, n Γ )F.
The method of the proof of [28, Theorem 1] is based on analysis of the closed system of differential equations  14)]). Here ν is the tangent angle such that t Γ = (cos ν, sin ν) T and n Γ = (sin ν, − cos ν) T . Since the arc-length parametrization can be expressed in terms of the relative local length r via ds = |∂ u X| du = rL(Γ t )( Φ(κ Γ ) /Φ(κ Γ ))du the system (21) is indeed closed.
The system of equations (21) can be rewritten as an abstract differential equation: The mapping H : E 1 → E 0 defined by the right-hand side of (21) is C 1 smooth on some open neighborhood of the initial condition Y ini ∈ E 1 . More precisely, Here the principal part H 0 (Y ) has the form: Indeed, sinceβ = aκ Γ − b − cF where a, b, c are given as in (12), (13), (14), we have Therefore all the terms ∂ s a, ∂ 2 s a, ∂ 2 s b, ∂ 2 s c can be expressed in terms of κ Γ , ∂ s κ Γ , the tangent angle ν and the position vector X.
Following the proof of [28, Theorem 1], the linearization H 0 (Ȳ ) of the principal part atȲ = (κ Γ ,ν,r, X) can be written in the form H 0 (Ȳ ) =D(Ȳ )∂ 2 u +C(Ȳ ), whereD (Ȳ ) = diag(āḡ −2 ,āḡ −2 , 0,āḡ −2 ) ∈ c ε (S 1 ) andC(Ȳ ) is a 5 × 5 matrix with coefficients smoothly depending onκ Γ ,ν,r, X and thus belonging to c ε (S 1 ). The rest of the proof is similar to that of [28,Theorem 1]. The principal partD∂ 2 u is a generator of an analytic semigroup on E 0 with the domain E 1 and it belongs to the maximal regularity class M(E 1 , E 0 ) on the pair of spaces (E 1 , E 0 ). Since H 1 ∈ C 1 (E 1 2 , E 0 ) the remaining operatorC + H 1 (Ȳ ) has the zero relative norm with respect to the principal part (c.f. [2]). Indeed, for the linear operator B =C + H 1 (Ȳ ) we have B ∈ L(E 1 2 , E 0 ). Taking into account the interpolation inequality between c ε and c 2+ε spaces we conclude the following inequality: E0 , where C 0 > 0 is a generic positive constant. Using Young's inequality ab ≤ δa 2 + b 2 /(4δ), δ > 0, we can conclude that the linear operator B (now considered as a linear operator from E 1 into E 0 ) has the relative zero norm, i.e. for any δ > 0 there exists a constant K δ > 0 such that BY E0 ≤ δ Y E1 + K δ Y E0 . With regard to [2, Lemma 2.5] the maximal regularity class M(E 1 , E 0 ) is closed with respect to perturbations by linear operators with the zero relative norm. Thus H (Ȳ ) ∈ M(E 1 , E 0 ). Applying the abstract theory due to Angenent [2, Theorem 2.7] the proof of the local existence and continuation follows. Remark 1. In Proposition 1 we assumed that the initial curve is C 4+ε smooth which might be considered as restrictive when compared to standard assumptions requiring less smoothness on initial data (c.f [12]). However, as the nontrivial curvature adjusted tangential velocity α depends on the curvature and tangent angle we had to consider the full governing system equations (21) including equations for the curvature, tangent angle. Following the methods based on maximal regularity we have to assume higher smoothness on the initial data as both the curvature and tangent angle should belong to the space c 2+ε (S 1 ) and so the initial curve is assumed to be C 4+ε .
Integrating equation (20) with respect to s, the tangential velocity α is a smooth mapping with values in c ε (S 1 ) provided that the curvature k belongs to the space c 1+ε (S 1 ) only. Following the theory due to Lunardi [18] (see also [2,Theorem 2.11]), one can expect that the optimal regularity assumption that the initial curve is C 3+ε smooth only. 6. Preservation of the surface area. The area A(Γ t ) enclosed by a curve Γ t evolved in the normal direction with the velocity v Γ satisfies the following identity: (c.f. [14], see also [16]). If the normal velocity is given by (3) we obtain d dt A(Γ t ) = 0, i.e. the enclosed area of planar curves Γ t , t ≥ 0, is preserved (c.f. Gage [14]). Let us consider evolution of closed curves G t , t ≥ 0, on a given surface M ⊂ R 3 according to the geometric equation (1). The surface M is prescribed as a graph of a function ϕ = ϕ(x), x ∈ Ω ⊂ R 2 , and the curve G t is parametrized by means of the vector valued mapping X(u, t) as G t = {(X, ϕ(X)) T : X ∈ Γ t }, where Γ t is the vertical projection of G t to R 2 . The resulting governing equations for the unknown parametrization are given by system (15) and with coefficients a, b, c given by (12), (13), (14).
Next we prove an analogous identity to (22) for evolving surface curves G t , t ≥ 0. As a consequence we will prove that the flow (1) preserves the surface area enclosed by G t . This area is denoted as A(G t ) and it can be expressed as follows: We extend the identity (1) to the case of surface curves G t evolving on the surface M given by function ϕ. Although the proof is straightforward, we present it for reader's convenience.
Proposition 2. Suppose {G t } t≥0 is a family of C 1 smooth, non-selfintersecting closed curves evolving on a surface M ⊂ R 3 with the normal velocity V G , which is a graph of a smooth function ϕ = ϕ(x). Then the identity holds for all t ≥ 0. In particular, if the normal velocity V G is given by (1) then the area A(G t ) = A(G ini ) is preserved for all t ≥ 0.
Proof. According to (23), the area of the surface on M enclosed by G t is given by A(G t ) = Int(Gt) dX = Int(Γt) 1 + |∇ϕ(x)| 2 dx. Let F 1 , F 2 : Ω → R be any smooth functions and such that One can easily construct such functions, e.g. F 1 = F 1 (x 1 ) and F 2 = F 2 (x 1 , x 2 ) is a primitive function to 1 + |∇ϕ(x 1 , x 2 )| 2 in the x 1 variable. Having such a pair of functions F = (F 1 , F 2 ) T , the area A(G t ) is given by Differentiating (25) with respect to t and taking into account that G t is a closed curve we obtain The terms containing the tangential velocity α canceled out because ∂ u X du = t Γ |∂ u X| du = t Γ ds.
Let us denote f ij = ∂Fi ∂xj , i, j = 1, 2. Then we calculate the scalar products in (26) as follows: Since dS = 1 + (∇ϕ · t Γ ) 2 ds and due to the fact that the normal velocity V G is related with β by (10), equation (26) becomes which is an analogy to the identity (22). Hence provided that V G = −K G + F, i.e. V G is given by the geometric equation (1). The surface enclosed by the curve G t remains constant for all t ≥ 0.
Proposition 3. The flow of surface curves G t , t ≥ 0, with the normal velocity V G given by (1) is the length shortening flow, i.e. d dt L(G t ) < 0 unless G t = G ini is a stationary curve having a constant geodesic curvature.
Proof. If the family of the surface curves G t , t ≥ 0, evolves in the outer normal direction by the velocity V G then, by [22,Eq. (32)], the length L(G t ) satisfies the identity: Clearly, the equality occurs if and only if K G is constant on G which implies V G = 0, i.e. G is a stationary curve.
7. Numerical approximation. The spatial discretization of (15) is based on the method of flowing finite volumes, which was applied and analyzed by Mikula anď Sevčovič in [21]. The principle of the method is that the discrete nodes X i = X(u i , t) for i = 0, . . . , M , X 0 = X M , and X 1 = X M +1 , and the dual nodes X i± 1 2 = X(u i± 1 2 , t) for u i± 1 2 = u i ± h/2, i = 1, . . . , M − 1, h = 1/M , are placed along the curve Γ t as illustrated in Figure 2.
Then parametric equations (15) are integrated over the dual segment between u i− 1 2 and u i+ 1 2 , which surrounds the discrete node X i , i = 1, . . . , M − 1. This integration results into the following expressions: Evaluation of the first integral on the right-hand side yields: where we have assumed a is constant on [u i− 1 2 , u i+ 1 2 ]. We denote the following discrete quantities where X 0 = X M and X M +1 = X 1 in the case of closed curve Γ t . Then, the curvature κ Γ is approximated as The approximation of the unit tangent a normal vectors is as follows: Finally, the discrete geodesic curvature (7) can be computed as: The quantities (12), (13), (14) from the parametric equation (15) are given by: and the nonlocal scalar force F becomes The area A(G t ) is calculated by means of equation (23) and the Green formula (see (25)). Having a pair of functions F 1 and F 2 such that The differences between the quantities A(G ini ) and A(G Ti ) are evaluated by means of the maximum norm both of them depending on the number M of finite volumes. Assuming that both error estimates are depending on the number of finite volumes as follows: then, the value of the Experimental order of convergence (EOC) between two levels of meshes with M 1 and M 2 finite volumes is given by: .
In Table 1 we summarize the parametric forms of initial conditions X ini and functions ϕ describing surfaces M for the following computational examples. The motion of the initial curve is driven by the normal velocity V G given by (1) and for each of the following examples we present how the curve G t asymptotically approaches the stable position, while the area enclosed by G t is preserved. The computations were performed for the number of segments in the flowing finite volume method ranging from 100 to 500. For each example we calculated EOC's which are summarized in their respective tables.
Example 1. The first example is presented in Figure 3. It illustrates, how an initial 5-leaf shape curve projected on the surface of the sphere evolves in the time interval   Table 1 is given as r(u) = (1 + 0.65 cos(10πu)). In Table 2 we present the values of EOC for the enclosed area.   Table 3 we present the values of EOC for the enclosed area.   Example 3. In the third example presented in Figure 5, evolution of the curve G t projected to the surface with sinus profile is illustrated. In this case, the time interval was [0,8]. With M = 200 finite volumes, we calculated that the initial curve G ini encloses the area of A(G ini ) = 7.15954, while the area enclosed at the final time T = 8 is A(G T ) = 7.15838. In Table 4 we present the values of EOC for the enclosed area.
Example 4. The fourth example is presented in Figure 6, and it shows evolution of the curve G t projected to the surface with saddle profile. In this case, the time interval was [0, 15]. With M = 200 finite volumes, we calculated that the initial curve G ini encloses the area of A(G ini ) = 2.30099, while the area enclosed at the final time T = 15 is A(G T ) = 2.30145. In Table 5 we present the values of EOC for the enclosed area.
9. Conclusion. We analyzed a non-local geometric flow preserving surface area enclosed by a closed curve on a given surface evolved in the normal direction by the geodesic curvature and the external force. We derived the form of the normal velocity of a nonlocal geometric flow that preserves the initial enclosed area. It was shown that the surface area preserving flow decreases the length of evolved  surface curves. Local existence and continuation of classical smooth solutions to the governing system of partial differential equations were also shown. As a numerical approximation scheme we proposed a method of flowing finite volumes method for spatial discretization in combination with the Runge-Kutta method for solving resulting system of ODEs. The scheme exhibited 2 nd order of experimental convergence. Several computational examples of evolution of surface curves were presented.