Variational discretization of axisymmetric curvature flows

We present natural axisymmetric variants of schemes for curvature flows introduced earlier by the present authors and analyze them in detail. Although numerical methods for geometric flows have been used frequently in axisymmetric settings, numerical analysis results so far are rare. In this paper, we present stability, equidistribution, existence and uniqueness results for the introduced approximations. Numerical computations show that these schemes are very efficient in computing numerical solutions of geometric flows as only a spatially one-dimensional problem has to be solved. The good mesh properties of the schemes also allow them to compute in very complex axisymmetric geometries.


Introduction
Numerical approximations of curvature flows such as the mean curvature flow and the Gauss curvature flow have been studied intensively during the last 30 years. In many situations the axisymmetry of these geometric flows can be used to reduce the dimension of the governing equations, and so numerical methods have been used frequently in such axisymmetric settings. However, results on the numerical analysis of such schemes so far are rare. In this paper we present parametric finite element approximations for axisymmetric curvature flows, and carefully analyse their properties.
In general, in curvature driven evolution equations the normal velocity of a hypersurface in R 3 is given by an expression involving the mean and/or the Gauss curvature of the surface. Evolving surfaces are of interest in geometry, and they can appear in application areas such as materials science, for example as grain boundaries. In addition, evolution laws involving the curvature of the surface arise in situations, where surface quantities are coupled to the surrounding volume by additional fields, which for example arises in the evolution of phase boundaries or in two-phase flow. In any case solving the evolution law for the surface with a stable discretization of curvature is a corner stone of a reliable and efficient numerical method.
Approaches to solve surface evolution equations numerically involve different descriptions of the evolving surface. Traditionally level set methods, phase field methods or parametric front tracking methods have been used and we refer to [17], and the references therein, for further information on numerical methods for geometric flows.
In this paper we aim to numerically compute a family of hypersurfaces (S(t)) t≥0 ⊂ R 3 , which we later assume to be axisymmetric, and which fulfills a geometric evolution law involving its principal curvatures. We will focus on the mean curvature flow, which for S(t) is given by the evolution law and which is the L 2 -gradient flow for H 2 (S(t)), since Here V S denotes the normal velocity of S(t) in the direction of the normal ν S(t) . Moreover, k m is the mean curvature of S(t), i.e. the sum of the principal curvatures of S(t), see [34] for an introduction to the mean curvature flow.
We also consider the nonlinear mean curvature flow where f : (a, b) → R with −∞ ≤ a < b ≤ ∞, is a strictly monotonically increasing continuous function, as well as the volume preserving variant Possible choices for f are f (r) = |r| β−1 r, β ∈ R >0 , (1.4a) or f (r) = −r −1 (1.4b) for the inverse mean curvature flow. These two choices have applications for example in image processing or in general relativity, see [39,31] and the references therein. Of course, (1.2) with f (r) = r (1.4c) collapses to (1.1).
If Ω(t) denotes the region enclosed by S(t), i.e. S(t) = ∂Ω(t), then the flow (1.3) is such that d dt L 3 (Ω(t)) = S(t) where here we assume that ν S is the outer normal to Ω(t) on S(t). This justifies the expression volume preserving flow. These flows are of interest in geometry and we refer to [29,4,13,27] for more information.
More generally, we can also consider flows of the form where k g = k 1 k 2 denotes the Gauss curvature of S(t), with k 1 and k 2 the two principal curvatures. Of course, (1.6) with F (r, s) = f (r) reduces to (1.2). On the other hand, the choice F (r, s) = −s leads to the Gauss curvature flow see e.g. [8, (1.14)], where in (1.7) we again assume that ν S is the outer normal to Ω(t) on S(t). Such flows have found considerable interest in geometry recently and we refer to [25,38,37,40] for more information. One reason why the Gauss curvature flow is of particular interest, is because this flow allows to study the fate of the rolling stones, see [2].
In this paper, we consider the case that S(t) is an axisymmetric surface, that is rotationally symmetric with respect to the x 2 -axis. We further assume that S(t) is made up of a single connected component, with or without boundary. Clearly, in the latter case the boundary ∂S(t) of S(t) consists of either one or two circles that each lie within a hyperplane that is parallel to the x 1 − x 3 -plane. For the evolving family of surfaces we allow for the following types of boundary conditions. A boundary circle may assumed to be fixed, it may be allowed to move vertically along the boundary of a fixed infinite cylinder that is aligned with the axis of rotation, or it may be allowed to expand and shrink within a hyperplane that is parallel to the x 1 − x 3 -plane. Depending on the postulated free energy, certain angle conditions will arise where S(t) meets the external boundary. If the free energy is just surface area, H 2 (S(t)), then a 90 • degree contact angle condition arises. We refer to Section 2 below for further details, in particular with regard to more general contact angles.
The dimensionally reduced formulation has several severe advantages both analytically as well as numerically. In analysis it has been used for example to study the onset of singularities, see [30,21,35,33,37] and other singularity formation mechanisms, see [11,12]. Numerically it leads to equations which are far easier to solve and at the same time problems with the mesh topology do not occur. Therefore, axisymmetric settings have been frequently used for numerical computations of surface evolutions. For example, graph formulations for axisymmetric geometric evolution laws have been considered in [14,16,18], while a finite difference approximation of a parametric description for the evolution of general axisymmetric surfaces has been studied in [36]. Hence the latter is closely related to the presented work, although we stress that it does not contain any numerical analysis. Moreover, also more complex problems such as for example two phase flows or biomembranes, in which also curvature effects play a role, have been treated in an axially symmetric setting. We refer to [24,41,28,15,43], and we expect that our approach will have an impact on such more complex evolutions as well. In terms of the numerical analysis for the approximation of axisymmetric surface evolutions only very few results have appeared in the literature so far, see e.g. [16,18] in the context of a graph formulation for the higher order curvature flows surface diffusion and Willmore flow, respectively. To the best of our knowledge, our paper contains the first stability results for fully discrete approximations of axisymmetric mean curvature flow. In addition, we will consider the numerical analysis of approximations for axisymmetric higher order flows, such as surface diffusion and Willmore flow, in the forthcoming article [10].
The present authors in the last ten years introduced parametric finite element methods for geometric evolution equations which have the property that the mesh generically behaves well during the evolution. We also refer to the recent work [22] for a method which also leads to good meshes. This is an advantage compared to earlier front tracking approaches in which often the meshes degenerated during the evolution such that the computations had to be stopped. In a series of papers, [5,6,7,8,9], we were able to analyze mesh properties and showed stability results. In particular, in two dimensions a semi-discrete version of the method led to equidistribution of mesh points. In this paper we introduce a parametric finite element method for the axisymmetric formulations of the surface evolution equations discussed above relying on ideas of our earlier work. However, a lot of new techniques have to be introduced stemming partly from the fact that close to the axis of rotation the equations, depending on the formulation, become either singular or degenerate, and partly because one has to decide how to deal with the equidistribution property. We will discuss several ways to handle these issues and will show stability, equidistribution, existence and uniqueness results for the new schemes. This paper is organised as follows. In Section 2 we introduce several weak formulations which will be crucial for the parametric finite element approximations introduced later. In Section 3 we derive semidiscrete, i.e. continuous in time discrete in space discretizations, and discuss stability and equidistribution properties. Section 4 is devoted to fully discrete schemes for which existence results are shown for linear as well as nonlinear variants as well as uniqueness results for linear schemes. In addition, we show stability for a fully discrete, mildly nonlinear discretization. Finally, we present several numerical results demonstrating that the majority of the schemes led to efficient, reliable results for mean curvature flow as well as for fully nonlinear curvature flows including its mass preserving variants.

Weak formulations
Let R/Z be the periodic interval [0, 1], and set We consider the axisymmetric situation, where x(t) : I → R 2 is a parameterization of Γ(t). Throughout Γ(t) represents the generating curve of a surface S(t) that is axisymmetric with respect to the x 2 -axis, see Figure 1. In particular, on defining Π 3 3 (r, z, θ) = (r cos θ, z, r sin θ) T for r ∈ R ≥0 , z ∈ R , θ ∈ [0, 2 π] and Π 3 2 (r, z) = { Π 3 3 (r, z, θ) : θ ∈ [0, 2 π)} , we have that Here we allow Γ(t) to be either a closed curve, parameterized over R/Z, which corresponds to S(t) being a genus-1 surface without boundary. In particular, we always assume that, for all t ∈ [0, T ], x(ρ, t) . e 1 > 0 ∀ ρ ∈ I \ ∂ 0 I , (2.2a) x(ρ, t) . e 1 = 0 ∀ ρ ∈ ∂ 0 I , is a disjoint partitioning of ∂I, with ∂ 0 I denoting the subset of boundary points of I that correspond to endpoints of Γ(t) attached to the x 2 -axis. Moreover, ∂ D I ∪ 2 i=1 ∂ i I denotes the subset of boundary points of I that model components of the boundary of S(t). Here endpoints in ∂ D I correspond to fixed boundary circles of S(t), that lie within a hyperplane parallel to the x 1 − x 3 -plane R × {0} × R. Endpoints in ∂ 1 I correspond to boundary circles of S(t) that can move freely along the boundary of an infinite cylinder that is aligned with the axis of rotation. Endpoints in ∂ 2 I correspond to boundary circles of S(t) that can expand/shrink freely within a hyperplane parallel to the we introduce the arclength s of the curve, i.e. ∂ s = | x ρ | −1 ∂ ρ , and set where (·) ⊥ denotes a clockwise rotation by π 2 . On recalling (2.1), we observe that the normal ν S on S(t) is given by ( ν(ρ, t) . e 1 ) cos θ ν(ρ, t) . e 2 ( ν(ρ, t) . e 1 ) sin θ    for ρ ∈ I , t ∈ [0, T ] , θ ∈ [0, 2 π) . (2.5) Similarly, the normal velocity V S of S(t) in the direction ν S is given by For the curvature κ of Γ(t) it holds that An important role in this paper is played by the surface area of the surface S(t), which is equal to (2.8) Often the surface area, A( x(t)), will play the role of the free energy in our paper. But for an open surface S(t), with boundary ∂S(t), we consider contact energy contributions which are discussed in [23], see also [9, (2.21)]. In the axisymmetric setting the relevant energy is given by ∂S ( x(p, t) . e 1 ) x(p, t) . e 2 + π ∂S , for p ∈ ∂ 1 I, denotes the change in contact energy density in the direction of − e 2 , that the two phases separated by the interface S(t) have with the infinite cylinder at the boundary circle of S(t) represented by x(p, t). Similarly, ∂S , for p ∈ ∂ 2 I, denotes the change in contact energy density in the direction of − e 1 , that the two phases separated by the interface S(t) have with the hyperplane R × {0} × R at the boundary circle of S(t) represented by x(p, t). These changes in contact energy lead to the contact angle conditions for all t ∈ (0, T ]. In most cases, the contact energies are assumed to be the same, so that ∂S = 0, which leads to 90 • contact angle conditions in (2.10a,b), and means that (2.9) collapses to (2.8). See [9] for more details on contact angles and contact energies. We note that a necessary condition to admit a solution to (2.10a) or to (2.10b) is that In addition, we observe that the energy (2.9) is not bounded from below if (p) For later use we note that ∂S [( x t (p, t) . e 1 ) x(p, t) . e 2 + ( x(p, t) . e 1 ) x t (p, t) . e 2 ] + 2 π Moreover, we recall that expressions for the mean curvature and the Gauss curvature of S(t) are given by x . e 1 and K S = −κ ν . e 1 x . e 1 on I , (2.13) respectively; see e.g. [15, (6)]. More precisely, if k m and k g denote the mean and Gauss curvatures of S(t), then k m = κ S (ρ, t) and k g = K S (ρ, t) on Π 3 2 ( x(ρ, t)) ⊂ S(t) , ∀ ρ ∈ I , t ∈ [0, T ] .
In the literature, the two terms making up κ S in (2.13) are often referred to as in-plane and azimuthal curvatures, respectively, with their sum being equal to the mean curvature. We note that combining (2.13) and (2.7) yields that x . e 1 ν . (2.14) A weak formulation of (2.14) will form the basis of our stable approximations for mean curvature flow and surface diffusion. Clearly, for a smooth surface with bounded mean curvature it follows from (2.14) that which, on recalling (2.4), is clearly equivalent to A precise derivation of (2.16) in the context of a weak formulation of (2.14) will be given in the Appendix A.
We note that (2.18b) weakly imposes (2.16) and (2.10a,b). We observe that (2.17a) degenerates for x . e 1 = 0, i.e. when ρ ∈ ∂ 0 I. Hence this degeneracy is balanced by the condition (2.15). In fact, on recalling (2.7) it holds that We remark that the weak formulation (A) is close in spirit to the weak formulations introduced in [6,7] for mean curvature flow. In particular, the tangential component of x t is not prescribed, which on the discrete level leads to an equidistribution property.
Similarly to (2.18b), we observe that (2.22b) weakly imposes (2.16) and (2.10a,b). We remark that the weak formulation (B) in some sense is close in spirit to the weak formulations introduced in [19,20] for mean curvature flow. In particular, the tangential component of x t is fixed to be zero, as the right hand side of (2.22a) is normal, recall (2.7).
Choosing χ = η = ( x . e 1 ) x t ∈ V ∂ in (2.22a,b), we obtain, similarly to (2.20), that We remark that it does not appear possible to mimic either (2.20) for (A) or (2.23) for (B) on the discrete level. Hence, in order to develop stable approximations, we investigate alternative formulations based on (2.14).
We note that the variational formulation for κ S in (2.25b) has previously been employed in [24, p. 124].
Choosing η = x t ∈ V ∂ in (2.24b) and χ = κ S in (2.24a), we obtain for the formulation (C), on recalling (2.12), that Similarly, choosing η = x t ∈ V ∂ in (2.25b) and χ = κ S in (2.25a), we obtain for the formulation (D) that We observe that (2.24b) and (2.25b) weakly impose (2.10a,b). But, in contrast to (2.18b) and (2.22b), it is not obvious that they also weakly impose (2.16), due to the presence of the degenerate weight x . e 1 . However, we show in the Appendix A that in fact they also weakly impose (2.16).
We also note that the formulation (C) is loosely related to (A), in the sense that the tangential component of x t is not prescribed. But in contrast to (A), discretizations of (C) cannot be shown to have an equidistribution property. In a similar way, the formulation (D) is loosely related to (B). However, we observe that the velocity x t is no longer purely in the normal direction. Finally, we observe that the variable κ S can be eliminated from (C), by choosing χ = ν . η in (2.24a) for η ∈ V ∂ , and then combining (2.24a) and (2.24b). Similarly, κ S can be eliminated from (D) by choosing χ = η in (2.25a) for η ∈ V ∂ , and then combining (2.25a) and (2.25b). We remark that the formulation (C), with the variable κ S , as well as the formulation (A), are useful with a view towards introducing numerical approximations of higher order flows, such as surface diffusion, see [10].

Nonlinear mean curvature flow
It is a simple matter to extend the formulations (A) and (C) to the nonlinear flow (1.2). In principle this can also be achieved for (B) and (D), but as the mean curvature needs to be recovered from the mean curvature vector, the resulting formulations are less natural. Hence we concentrate on (A) and (C). For the former, replacing the right hand side in x . e 1 ) χ | x ρ | dρ yields a weak formulation for (1.2), which we call Similarly to (2.26), and using the same choices of η and χ, it can be shown that solutions to (C f ) satisfy which yields stability if f is monotonically increasing with f (0) = 0.
Finally, we may also generalize these nonlinear formulations to the volume preserving flow (1.3).
Choosing χ = 2 π in (2.31a) yields (2.30), as before. Moreover, and similarly to (2.28), it can be shown for solutions of (C f,V ) in the case (1.4c) that (2.32) where we have noted the Cauchy-Schwarz inequality. It does not appear possible to extend the stability result (2.32) to the case of more general f .

Semidiscrete schemes
, be a decomposition of [0, 1] into intervals given by the nodes q j , I j = [q j−1 , q j ]. For simplicity, and without loss of generality, we assume that the subintervals form an equipartitioning of [0, 1], i.e. that Clearly, if I = R/Z we identify 0 = q 0 = q J = 1.

The necessary finite element spaces are given by
For later use, we let π h : C(I) → V h be the standard interpolation operator at the nodes {q j } J j=0 . Let (·, ·) denote the L 2 -inner product on I, and define the mass lumped L 2 -inner product (f, g) h , for two piecewise continuous functions, with possible jumps at the nodes {q j } J j=1 , via 2) naturally extends to vector valued functions. It is easily shown that Assuming that | X h ρ | > 0 almost everywhere on I, and similarly to (2.4), we set We note that For later use, we let ω h ∈ V h be the mass-lumped Recall from (2.8) and (2.9) that We have, similarly to (2.12), that (3.8)

Mean curvature flow
In view of the degeneracy on the right hand side of (2.17a), and on recalling (2.19) and Our semidiscrete finite element approximation of (A), (2.18a,b), is given as follows.
The equidistribution property (3.11) can be shown by choosing We also remark that it follows from (3.6) that We note that mass lumping in (3.10b) is crucial for the proof of the equidistribution property (3.11). Hence we only consider the variant (A h ) h with mass lumping. Of course, in the case ∂ 0 I = ∅, an alternative scheme to (3.10a,b) is together with (3.10b). Note that if ∂ 0 I = ∅ then (3.10a) collapses to (3.13) with ν h replaced by ω h . Unfortunately, neither choice appears to lead to a stability proof.
In an attempt to prove stability, Moreover, considering for simplicity the case ∂ 0 I = ∅, and choosing Unfortunately, the right hand sides in (3.14) and (3.15) are not equal, recall (3.6), and so combining (3.14) and (3.15) does not yield a stability result. On the other hand, the function ( is not well defined, and cannot be chosen as a test function in (3.10a) or (3.13).
However, the fully discrete variant of (A h ) h , (3.10a,b), performs very well in practice.
The rescaling factor | ω h (q j , t)| 2 in (3.17) normalizes the discrete vertex normals ω h (q j , t), recall (3.12), which is the most natural approach. Similarly to (A h ) h , it does not appear possible to prove a stability result for (B h ) h .
However, it turns out that approximations of the formulations (C) and (D) can be shown to be stable. In particular, our semidiscrete approximations of (C), (2.24a,b), and (D), (2.25a,b), are given as follows, where we first define Here and throughout we use the notation · (h) to denote an expression with or without the superscript h, and similarly for the subscripts · ∂ 0 . I.e. the scheme (C h ) h employs mass lumping, recall (3.2), and seeks κ S (t) ∈ W h ∂ 0 , while the scheme (C h ) employs true integration throughout and seeks κ S (t) We observe that (C h ) h and (D h ) h do not depend on the values of κ h S and κ h S , respectively, on ∂ 0 I. Hence we fix these values to be zero by requiring that κ h S ∈ W h ∂ 0 and κ h S ∈ W h ∂ 0 , and by using a reduced set of test functions in (3.18a) and (3.19a). As a consequence, it scheme stability implicit tangential motion equidistribution yes no no Table 1: Properties of the different semidiscrete schemes for mean curvature flow.
seems at first that X h t is not defined on ∂ 0 I. However, X h on ∂ 0 I is determined through (3.18b) and (3.19b), respectively.
We have on choosing (3.19a) and η = X h t in (3.18b) and (3.19b), on recalling (3.8), that respectively, where we have recalled (3.6). This shows that they can be interpreted as natural L 2 -gradient flows of (3.7).
We observe that it is possible to eliminate κ h S from the schemes (D h ) (h) , which yields (3.19b) with κ h S replaced by X h t . Similarly, κ h S can be removed from the scheme (C h ) h to yield (3.18b) with κ h S ν h replaced by ( X h t . ω h ) ω h , on recalling (3.6). For the scheme (C h ) this elimination procedure is not possible.
For the reader's convenience, Table 1 summarises the main properties of all the schemes introduced in this section. given by (3.10a,b) with the right hand side in (3.10a) replaced by

Nonlinear mean curvature flow
These two schemes inherit the equidistribution property, recall (3.11). Replacing κ h S with π h [f (κ h S )] in (3.18a) yields the schemes (C f h ) (h) and similarly we can define (C f,V h ) (h) by replacing the right hand side in (3.18a) by .
(3.22) Similarly to (3.20), and using the same choices of η and χ, it can be shown that solutions to , which yields a stability bound for (C f h ) h if f is monotonically increasing with f (0) = 0. Of course, (3.22) is a discrete analogue of (2.28). Moreover, solutions to (C f,V h ) (h) , in the case (1.4c), satisfy similarly to (2.32), where here we have also used a Cauchy-Schwarz inequality for the mass lumped inner product (3.2). Finally, solutions to the scheme (C f,V h ) conserve the volume of the domain Ω h (t) ⊂ R 3 that is enclosed by the three-dimensional axisymmetric surface S h (t) that is generated by the curve Γ h (t). To see this, choose χ = 2 π in (3.18a), with the modified right hand side (3.22), to obtain recall (1.5). Here V h S h (t) denotes the normal velocity of S h (t) in the direction of ν h S h (t), the outer normal to Ω h (t) on S h (t), where ν h S h (t) is induced by ν h through a discrete analogue of (2.5). Using the same testing procedure for the scheme (C f,V h ) h yields that and so the enclosed volume is only approximately preserved, compare with (3.23). Finally, choosing χ = X h . e 1 in (A f,V h ) h , recall (3.21), also yields (3.24), and so an approximate volume preservation property.

Fully discrete schemes
Let ω m ∈ V h be the natural fully discrete analogue of ω h ∈ V h , recall (3.6), i.e.

Mean curvature flow
Similarly to (3.9), and given a κ m+1 ∈ V h , we introduce K m (κ m+1 ) ∈ V h such that Then our fully discrete analogue of (A h ) h , (3.10a,b), is given as follows. (4.3b) We make the following mild assumptions.
Note that the assumption (B) h , on recalling (3.6), is equivalent to assuming that dim span{ ω m (q j )} J j=0 = 2.
Proof. We note that since X m ∈ V h ∂ 0 satisfies the assumption (A), the right hand side of (4.3a) is well-defined. As (4.3a,b) is linear, existence follows from uniqueness. To investigate the latter, we consider the system: where we recall from (4.2) that λ ∈ V h with λ(q j ) = 1 q j ∈ I \ ∂ 0 I , 2 q j ∈ ∂ 0 I . Choosing χ = κ ∈ V h in (4.4a) and η = δ X ∈ V h ∂ in (4.4b) yields that It follows from (4.6) that κ = 0 and that δ X ≡ X c ∈ R 2 ; and hence that It follows from (4.7) and assumption (B) h that X c = 0. Hence we have shown that We remark that a fully discrete approximation of (B h ) h , (3.16a,b), is given by: In practice the scheme (4.8a,b), for reasonable time step sizes, can lead to oscillations and poor results, see e.g. Figure 3 below.

Lemma. 4.2.
Let X m ∈ V h ∂ 0 satisfy the assumption (A). There exists a unique solution Proof. Similarly to the proof of Lemma 4.1, we obtain that where (δ X, κ) ∈ V h ∂ × V h solve the linear homogeneous system corresponding to (4.8a,b). It follows from (4.9) that κ = 0 and then from the homogeneous variant of (4.8a) that δ X = 0. Hence we have shown that (4.8a,b) has a unique solution (δ Our fully discrete analogues of the schemes (C h ) (h) , (3.18a,b), and (D h ) (h) , (3.19a,b), are given as follows.
For the second variant, which is going to lead to systems of nonlinear equations and for which a stability result can be shown, we introduce the notation [r] ± = ± max{±r, 0} for r ∈ R. . We note that in practice we solve the nonlinear schemes with a Newton iteration, which in all our experiments always converged with at most three iterations.

Remark. 4.3.
Similarly to the semidiscrete variants, we observe that in most of the above fully discrete schemes it is possible to eliminate the discrete curvatures, κ m+1 S or κ m+1 S . For example, on recalling (3.6) and on choosing χ = π h [ η . ω m ] ∈ W h (∂ 0 ) in (4.11a) for η ∈ V h ∂ , the scheme (C m, ) h reduces to: Find δ X m+1 ∈ V h ∂ such that, for all η ∈ V h ∂ , ∂S ( X m+1 (p) . e 1 ) η(p) . e 3−i . (4.14) and similarly for (C m ) h , (D m ) (h) and (D m, ) (h) , with the latter leading to (4.13b) with κ m+1 S replaced by (∆t m ) −1 ( X m+1 − X m ). For the schemes (C m ) and (C m, ) this elimination procedure is not possible. A related variant to (4.14) is given by: We make the following mild assumption.
Lemma. 4.4. Let X m ∈ V h ∂ 0 satisfy the assumptions (A) and (C (∂ 0 ) ) (h) . Then there exists a unique solution (δ X m+1 , κ m+1 Proof. As (4.10a,b) is linear, existence follows from uniqueness. To investigate the latter, we consider the system: It immediately follows from (4.17) and the assumption (A) that κ S = 0, and that δ X ≡ X c ∈ R 2 . Hence it follows from (4.16a) that X c . z = 0 for all z ∈ Z (h) (∂ 0 ) , and so assumption (C (∂ 0 ) ) (h) yields that X c = 0. Hence we have shown that (C m ) (h) has a unique solution (δ X m+1 , κ m+1 Lemma. 4.5. Let X m ∈ V h ∂ 0 satisfy the assumption (A). Then there exists a unique solution (δ X m+1 , κ m+1 Proof. Similarly to the proof of Lemma 4.4, we obtain that solve the linear homogeneous system corresponding to (4.12a,b). It immediately follows from (4.18) and the assumption (A) that κ S = 0, and that δ X = X c ∈ R 2 . Combined with the homogeneous variant of (4.12a) these imply that X c = 0. Hence we have shown that (4.12a,b) has a unique solution (δ X m+1 , κ m+1 Theorem. 4.6. Let ∂ 0 I = ∅ and let X m ∈ V h satisfy the assumption (A). Then there exists a solution (δ X m+1 , κ m+1 Upon eliminating κ m+1 S from (D m, ) (h) , we can rewrite it as: Given On recalling assumption (A), we note that ∂ 0 I = ∅ implies that µ = min I X m . e 1 > 0. It holds that ∂S | X m (p) . e 1 | η(p)| . (4.20) In relation to the second term on the right hand side of (4.20) we observe, on recalling (3.3), that Moreover, it holds that

22) recall (3.3), and similarly
where in the last inequality we have used a Young's inequality and observed from (3.2) and for p ∈ ∂ i I, i = 1, 2. Now choosing ∆t m < 3 µ 2 in (4.24) we obtain that holds for γ sufficiently large, and so the existence of a solution δ X then follows immediately from (4.13a).
We remark that although one can eliminate κ m+1 S from the scheme (C m, ) h , recall (4.14), one cannot adapt the above proof for the scheme (D m, ) h , as one would obtain (4.20) with (| η| 2 , | X m ρ |) h replaced by (| η . ω m | 2 , | X m ρ |) h , and so it is no longer possible to bound e.g. the second term on the right hand side of (4.20).
We note that the scheme (4.15) can also be shown to be unconditionally stable, i.e. a solution to (4.15) satisfies

Nonlinear mean curvature flow
It is a simple matter to extend the presented fully discrete approximations to the nonlinear mean curvature flow (1.2) and the volume preserving variant (1.3). We recall the fully 3d parametric finite element schemes (2.4a,b) and (2.5) in [7] for the approximation of (1.2) and (1.3), respectively. We can now define their natural axisymmetric analogues. For example, the natural adaptation (A f m ) h of the scheme (A m ) h to (1.2) is given by (4.3a,b) with (4.3a) replaced by Similarly, the natural adaptation (C f m, ) (h) of the scheme (C m, ) (h) to (1.2) is given by (4.11a,b) with (4.11a) replaced by . (4.28) Similarly to (4.25), with the same choices of χ and η, it is then possible to prove that This is linear scheme for which existence of a unique solution, provided that the assumption (B) h holds, can easily be shown. Moreover, solutions to the semidiscrete variant of (A F m ) h satisfy the equidistribution property (3.11).
We also consider the ratio between the longest and shortest element of Γ m , and are often interested in the evolution of this ratio over time.

5.1
Numerical results for mean curvature flow

Sphere
It is easy to show that a sphere of radius r(t), with is a solution to (1.1). We use this true solution for a convergence test for the various schemes for mean curvature flow, similarly to Table 1 in [7]. Here we start with a nonuniform partitioning of a semicircle of radius r(0) = r 0 = 1 and compute the flow until time T = 0.125. In particular, we have ∂ 0 I = ∂I = {0, 1} and we choose , j = 0, . . . , J ,      Table 1 in [7], we see that the errors for the axisymmetric schemes are significantly smaller than the values in Table 1 in [7] for similar discretization parameters. It is clear from Table 2 that the schemes (A m ) h and (B m ) h appear to converge with the optimal convergence rate of O(h 2 ). Similarly, Tables 3 and 4 suggest that the schemes (C m, ) (h) , (D m ) (h) and (D m, ) (h) converge with an order slightly less than quadratic. We note that the linear schemes (C m ) (h) lead to solutions X m with min ρ∈I X m (ρ) . e 1 < 0, and so we cannot complete the evolutions. In particular, in practice the two boundary elements shrink in size due to the scheme's tangential motion. Once the element has shrunk to a length almost zero, the freely moving vertex can become negative. The linear schemes (D m ) (h) behave well, on the other hand. But as they are very close to the nonlinear schemes (D m, ) (h) , from now on we concentrate on the schemes (D m, ) (h) , (C m, ) (h) , (A m ) h and (B m ) h .

Torus
We repeat the two torus experiments in Figures 5 and 6 in [7]. To this end, we let ∂I = ∅. For an initial torus with radii R = 1, r = 0.7, we obtain a surface that closes up towards    Figure 2 for the simulation results for the scheme (A m ) h , for the discretization parameters J = 256 and ∆t = 10 −4 . On the other hand, for an initial torus with radii R = 1, r = 0.5, we obtain a shrinking evolution towards a circle, as in [7, Fig. 6]. See Figure 3 for the evolution for the scheme (A m ) h , again for J = 256 and ∆t = 10 −4 . On repeating the numerical experiment for the (B m ) h , we observe strong oscillations, as shown on the right of Figure 3. These oscillations become smaller in magnitude as ∆t is decreased. The remaining schemes can integrate the evolution shown in Figure 3 in a stable way, and their numerical results are very close to the ones displayed in Figure 3 for the scheme (A m ) h . However, the schemes differ in the exhibited tangential motions, which leads to diverse evolutions of the ratio r m , see Figure 4. The best distribution of mesh points is shown by the schemes (A m ) h and (C m, ) h , followed by (C m, ). The most nonuniform distribution of mesh points can be observed for the two schemes (D m, ) (h) .

Cylinder
For the scheme (A m ) h we repeat the singular evolution from [21, Fig. 1]. To this end, we set ∂ D I = ∂I = {0, 1}. In particular, starting with a cylinder, mean curvature flow leads to a pinch-off. We show the results for the scheme (A m ) h , with the discretization parameters J = 128 and ∆t = 10 −4 , in Figure 5.
For the next two experiments, we consider a cylinder attached to two parallel hyperplanes, with prescribed contact angle conditions, recall (2.10b). To this end, we set ∂ 2 I = ∂I = {0, 1} and use the discretization parameters J = 128 and ∆t = 10 −3 . Letting and starting with a cylinder, the evolution yields a growing catenoid-like surface, see Figure 6. We observe convergence to a travelling wave type solution, with the associated energy unbounded from below. We conjecture that the profile of the curve approaches in the limit the so-called grim reaper solution, see [34, p. 15] and [26], g(ρ, t) = (z 0 + π 3 t) e 1 + (− 3 π ln cos π 3 (ρ − 1 2 ) , ρ) T , (5.5)  where z 0 ∈ R specifies the position of the travelling wave solution at time t = 0. In fact, plotting g(ρ, t) − (z 0 + π 3 t) e 1 at time t = 4 versus X m − (min ρ∈I X m . e 1 ) e 1 for our final solution in Figure 6, yields perfect agreement between the two graphs. We conjecture that the speed of the travelling wave type solution will approach π 3 , the speed of (5.5). To test this conjecture, we continue the evolution until t = 100 and plot the evolution of X m (0) . e 1 over time, comparing the graph with a suitably chosen line with slope π 3 , see Figure 7. As we can see, the speed of the curve does indeed approach π 3 .
If we let ∂S = 1 2 , on the other hand, we observe a shrinking surface, with the radius of the contact circles eventually converging to zero. On reaching two single contact points with the external substrates, we allow the discrete surface to detach from the two hyperplanes and to continue the evolution as a closed genus 0 surface, see Figure 8 for the evolution. To allow for an accurate resolution of the detaching, we employ the smaller

Surface patch within a cylinder
For the next experiment, we consider a disk attached to an infinite cylinder of radius 1, with prescribed contact angle conditions, recall (2.10a). To this end, we set ∂ 0 I = {0} and ∂ 1 I = {1}, and use the discretization parameters J = 128 and ∆t = 10 −3 . Letting (1) ∂S = − 1 2 and starting with a disk, the evolution seems to converge to a translating surface patch, see Figure 9. Taking the angle condition (2.10a), there is a unique convex scaled surface grim reaper profile moving with constant speed by translation. We conjecture that a general class of initial data will converge to this shape for large times. We refer to [1] for more information on the grim reaper analogues in higher dimensions.   Figure 10. An insight that we gain from this set of experiments is that the tangential motion displayed by the scheme (C m, ) h can lead to very nonuniform meshes. Hence, for the remainder of this paper, we will only present numerical results for the two schemes (A m ) h and (C m, ) and their nonlinear variants. Note that the former is a linear fully discrete approximation of (A h ) h , for which the equidistribution property (3.11) holds. The latter, on the other hand, is a nonlinear scheme that is unconditionally stable, recall Theorem 4.7. As the results for (A m ) h and (C m, ) are often indistinguishable, we only visualize the numerical results for the former from now on.

Genus 0 surface
An experiment for a cigar shape can be seen in Figure 11. Here we have once again that ∂ 0 I = ∂I = {0, 1}. The discretization parameters are J = 128 and ∆t = 10 −4 . The relative volume loss for this experiment for (A f,V m ) h is 0.09%, while for (C f,V m, ) it is −0.01%. The same experiment for the scheme (C f,V m, ) h yields a very nonuniform mesh, with the final ratio r m > 430. An experiment for a disc shape is shown in Figure 12. The discretization parameters are J = 128 and ∆t = 10 −4 . The relative volume loss for this experiment for (A f,V m ) h is −0.02%, while for (C f,V m, ) it is −0.01%. Once again, the scheme (C f,V m, ) h yields a very nonuniform mesh for this simulation, with the final ratio r m > 145.

Genus 1 surface
We repeat the simulation in Figure 3 for conserved mean curvature flow, i.e. (1.3) with (1.4c), using the scheme (A f,V m ) h . Conservation of the enclosed volume means that the torus can no longer shrink to a circle. Hence the torus now attempts to close up and change topology, as can be seen from the numerical results in Figure 13. As for the original experiment, we use the discretization parameters J = 256 and ∆t = 10 −4 . The relative enclose volume loss for this experiment is −0.00%. The evolutions for the schemes (C f,V m, ) h and (C f,V m, ) are nearly identical to what is shown in Figure 13, with a relative volume loss Figure 11: (A f,V m ) h for (1.4c). Conserved mean curvature flow for a cigar. Plots are at times t = 0, 0.1, . . . , 1. We also visualize the axisymmetric surface S m generated by Γ m at time t = 0.3. On the right are plots of the discrete energy and the ratio r m and, as a comparison, a plot of the ratio r m for the scheme (C f,V m, ).
of 0.01% in both cases.
Finally, we present an example for conserved mean curvature flow, (1.3) with (1.4c), for the scheme (A f,V m ) h with the initial data X 0 parameterizing a closed spiral, so that the approximated surface has genus 1. As can be seen from Figure 14, the spiral slowly untangles, until the surface becomes a torus. For this experiment we use the discretization parameters J = 1024 and ∆t = 10 −6 . The relative enclosed volume loss for this experiment is 0.01%. The evolutions for the schemes (C f,V m, ) h and (C f,V m, ) are nearly identical to what is shown in Figure 14, with a relative volume loss of 0.01% in both cases.

Numerical results for nonlinear mean curvature flow
Similarly to (5.2), it is easy to show that a sphere of radius r(t), with is a solution to (1.2) with (1.4a). We use this true solution for a convergence test for β = 1 2 , similarly to Table 2 Table 5.
We repeat the same convergence experiment for the inverse mean curvature flow, where we note that a sphere of radius r(t), with is a solution to (1.2) with (1.4b). The errors are reported in Table 6. We recall that these numbers can be compared to the fully 3d results in Table 3 in [7]. It is clear from Tables 5 and 6 that the solutions to the scheme (A f m ) h appear to converge with the optimal Figure 13: (A f,V m ) h for (1.4c). Conserved mean curvature flow for a torus with radii R = 1, r = 0.5. Plots are at times t = 0, 0.01, . . . , 0.14. We also show a plot at time t = 0.145, together with a plot of the discrete energy over time. Below we visualize the axisymmetric surface S m generated by Γ m at times t = 0 and t = 0.145.  .  , on the other hand, the solutions appear to converge with an order less than quadratic, and closer to 3 2 . We believe that these lower convergence rates are caused by the nonuniform meshes induced by the scheme (C f m, ), recall Figure 10, and also by the degeneracy of the coefficients x . e 1 in (C f ).
In the next experiment we repeat the simulation in [7, Fig. 8] for the inverse mean curvature of a torus with radii R = 1, r = 0.25. We recall that for this nonconvex initial data, with κ S (·, 0) < 0, the classical inverse mean curvature develops a singularity in finite time, see also [25,40]. For the axisymmetric setting we use I = R/Z, so that ∂I = ∅. As the discretization parameters for the scheme (A f m ) h we use J = 256 and ∆t = 10 −4 . See Figure 15 for the simulation results. Similarly to the results in [7, Fig.  8], the discrete solution becomes unphysical after around time 0.52, where we conjecture that the singularity for the continuous flow occurs. Figure 15: (A f m ) h for (1.4b). Inverse mean curvature flow for a torus with radii R = 1, r = 0.25. Plots are at times t = 0, 0.05, . . . , 0.55. We also visualize the axisymmetric surface S m generated by Γ m at times t = 0 and t = 0.5.

Numerical results for Gauss curvature flow
An experiment for Gauss curvature flow, (1.7), for the same initial data as in Figure 11, can be seen in Figure 16. Here we have once again that ∂ 0 I = ∂I = {0, 1}. The discretization parameters for the scheme (A F m ) h from Remark 4.8 are J = 128 and ∆t = 10 −5 . As a comparison, we also show the evolution for standard mean curvature flow, computed with the scheme (A m ) h , in Figure 16.
The Gauss curvature flow is not well-defined for general hypersurfaces as for nonconvex hypersurfaces the resulting equation is not parabolic, see [3,32]. In the axisymmetric situation the degrees of freedom are reduced and the resulting equation is x . e 1 on I , which is parabolic as long as ν . e 1 is positive. In conclusion, even in the axisymmetric case the evolution is not well-defined if the initial surface has the topology of a torus, and so we do not present results for genus 1 surfaces. However, although this would not be possible in the general formulation we can start the Gauss curvature flow in the axisymmetric case with some nonconvex initial data, and we do so in the simulation in Figure 17, where we used the discretization parameters J = 128 and ∆t = 10 −5 . For a mathematical analysis for Gauss curvature flow in the axisymmetric case we refer to [32].

Conclusions
We have derived and analysed various numerical schemes for the parametric approximation of axisymmetric mean curvature flow, its nonlinear and volume conserving variants, as well as more general curvature flows. The main fully discrete schemes to consider for standard mean curvature flow are (A m ) h and (C m, ). Here we have dismissed the scheme (B m ) h , as it can show oscillations in practice, recall Figure 3, as well as the scheme (C m, ) h , Figure 16: (A F m ) h for (1.7). Gauss curvature flow for a cigar, on the left. Plots are at times t = 0, 0.01, . . . , 0.12. We also visualize the axisymmetric surface S m generated by Γ m at time t = 0.12. As a comparison, we show the evolution of (A m ) h for mean curvature flow, (1.1), on the right. Here the plots are at times t = 0, 0.01, . . . , 0.06, and we visualize the axisymmetric surface S m generated by Γ m at time t = 0.06.  7). Gauss curvature flow for nonconvex initial data. Plots are at times t = 0, 0.05, . . . , 0.6. We also visualize the axisymmetric surface S m generated by Γ m at times t = 0, t = 0.2 and t = 0.6. as it can display very nonuniform meshes in simulations where the discrete curves are attached to the x 2 -axis, recall Figure 10 for its variant (C f,V m, ) h . We also do not consider the schemes (D m, ) (h) , as they have no advantage over (C m, ) and as they can also exhibit very nonuniform meshes. Of the two schemes we consider, the scheme (A m ) h is a linear scheme that asymptotically leads to an equidistribution of mesh points, recall Remark 3.1. In addition, even though there is no stability proof for (A m ) h , in practice the discrete energy is always monotonically decreasing. The scheme (C m, ), on the other hand, is a nonlinear scheme that is unconditionally stable. The nonlinearity is only very mild, and so a Newton solver never takes more than 3 iterations in practice. Moreover, the distribution of vertices for (C m, ) may be worse than for (A m ) h , but coalescence of vertices is not observed in practice. Similar statements hold for the nonlinear variants (A f m ) h , (A f,V m ) h , (C f m, ) and (C f,V m, ), where the two conserving schemes show very good volume conservation properties in practice. Finally, for general curvature flows of the form (1.6), we propose the linear scheme (A F m ) h , which asymptotically exhibits equidistributed mesh points.
It follows that (B.2) is the Euler-Lagrange variation of the minimization problem: where Φ ∈ C 1 (R) denotes an antiderivative of f −1 . We note that Φ : R → R is strictly convex with Φ (f (0)) = f −1 (f (0)) = 0 and hence we obtain that Φ is bounded from below and is coercive.
In the following we establish that the continuous functional J h : V h ∂ → R is coercive, i.e. that J h ( η) → ∞ as η → ∞, where · is a fixed norm on V h . The main task is to bound the growth of the linear term m ( η) in terms of the first two terms in (B.3b).