Abstract
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.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 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 \({{\mathbb {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. For example, parametric finite element approximations of curvature flows have been considered in [7, 20, 23, 35]. We refer to the review paper [18], and the references therein, for further information on numerical methods for general geometric evolution equations.
In this paper we aim to numerically compute a family of hypersurfaces \(({\mathscr {S}}(t))_{t \ge 0}\)\(\subset {{\mathbb {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 \({\mathscr {S}}(t)\) is given by the evolution law
and which is the \(L^2\)-gradient flow for \({\mathscr {H}}^2({\mathscr {S}}(t))\), since
for surfaces without boundary. Here \({\mathscr {V}}_{{\mathscr {S}}}\) denotes the normal velocity of \({\mathscr {S}}(t)\) in the direction of the normal \(\vec {\mathrm{n}}_{{\mathscr {S}}}\). Moreover, \(k_m\) is the mean curvature of \({\mathscr {S}}(t)\), i.e. the sum of the principal curvatures of \({\mathscr {S}}(t)\), see [37] for an introduction to the mean curvature flow.
We also consider the nonlinear mean curvature flow
where \(f{:}\,(a,b)\rightarrow {\mathbb {R}}\) with \(-\,\infty \le a<b\le \infty \), is a strictly monotonically increasing continuous function, as well as the volume preserving variant
Possible choices for f are
or
for the inverse mean curvature flow. These two choices have applications for example in image processing or in general relativity, see [33, 42] and the references therein. Of course, (1.2) with
collapses to (1.1).
If \(\varOmega (t)\) denotes the region enclosed by \({\mathscr {S}}(t)\), i.e. \({\mathscr {S}}(t) = \partial \varOmega (t)\), then the flow (1.3) is such that
where here we assume that \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) is the outer normal to \(\varOmega (t)\) on \({\mathscr {S}}(t)\). This justifies the expression volume preserving flow. These flows are of interest in geometry and we refer to [4, 14, 29, 31] 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 \({\mathscr {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\), for closed surfaces, leads to the Gauss curvature flow
see e.g. [8, (1.14)], where in (1.7) we again assume that \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) is the outer normal to \(\varOmega (t)\) on \({\mathscr {S}}(t)\). Such flows have found considerable interest in geometry recently and we refer to [27, 40, 41, 43] 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 \({\mathscr {S}}(t)\) is an axisymmetric surface, that is rotationally symmetric with respect to the \(x_2\)-axis. We further assume that \({\mathscr {S}}(t)\) is made up of a single connected component, with or without boundary. Clearly, in the latter case the boundary \(\partial {\mathscr {S}}(t)\) of \({\mathscr {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 \({\mathscr {S}}(t)\) meets the external boundary. If the free energy is just surface area, \({\mathscr {H}}^2({\mathscr {S}}(t))\), then a \(90^\circ \) degree contact angle condition arises. We refer to Sect. 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 [22, 32, 36, 38, 40] and other singularity formation mechanisms, see [12, 13]. 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 [15, 17, 19], while a finite difference approximation of a parametric description for the evolution of general axisymmetric surfaces has been studied in [39]. 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 [16, 26, 30, 44, 45], 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. [17, 19] 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 consider the numerical analysis of approximations for axisymmetric higher order flows, such as surface diffusion and Willmore flow, in the recently appeared article [11].
The present authors in the last 10 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 [23] 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 Sect. 2 we introduce several weak formulations which will be crucial for the parametric finite element approximations introduced later. In Sect. 3 we derive semidiscrete, i.e. continuous-in-time discrete in space discretizations, and discuss stability and equidistribution properties. Section 4 is devoted to the fully discrete schemes that turn out to be most practical. For the linear schemes we prove existence and uniqueness results, while a stability result can be shown for a mildly nonlinear discretization. Additional fully discrete variants, and some of their properties, are summarized in “Appendix C”. Finally, we present several numerical results demonstrating that the majority of the schemes lead to efficient, reliable results for mean curvature flow as well as for fully nonlinear curvature flows including its mass preserving variants.
2 Weak formulations
Let \({{\mathbb {R}}} \slash {{\mathbb {Z}}}\) be the periodic interval [0, 1], and set
We consider the axisymmetric situation, where \(\vec {x}(t){:}\,\overline{I} \rightarrow {{\mathbb {R}}}^2\) is a parameterization of \(\varGamma (t)\). Throughout \(\varGamma (t)\) represents the generating curve of a surface \({\mathscr {S}}(t)\) that is axisymmetric with respect to the \(x_2\)-axis, see Fig. 1. In particular, on defining
and
we have that
Here we allow \(\varGamma (t)\) to be either a closed curve, parameterized over \({{\mathbb {R}}} \slash {{\mathbb {Z}}}\), which corresponds to \({\mathscr {S}}(t)\) being a genus-1 surface without boundary. Or \(\varGamma (t)\) may be an open curve, parameterized over [0, 1]. Then \(\varGamma (t)\) has two endpoints, and each endpoint can either correspond to an interior point of \({\mathscr {S}}(t)\), or to a boundary circle of \({\mathscr {S}}(t)\). Endpoints of \(\varGamma (t)\) that correspond to an interior point of the surface \({\mathscr {S}}(t)\) are attached to the \(x_2\)-axis, on which they can freely move up and down. For example, if both endpoints of \(\varGamma (t)\) are attached to the \(x_2\)-axis, then \({\mathscr {S}}(t)\) is a genus-0 surface without boundary. If only one end of \(\varGamma (t)\) is attached to the \(x_2\)-axis, then \({\mathscr {S}}(t)\) is an open surface with boundary, where the boundary consists of a single connected component. If no endpoint of \(\varGamma (t)\) is attached to the \(x_2\)-axis, then \({\mathscr {S}}(t)\) is an open surface with boundary, where the boundary consists of two connected components.
In particular, we always assume that, for all \(t \in [0,T]\),
where \(\partial _D I \cup \bigcup _{i=0}^2 \partial _i I = \partial I\) is a disjoint partitioning of \(\partial I\), with \(\partial _0 I\) denoting the subset of boundary points of I that correspond to endpoints of \(\varGamma (t)\) attached to the \(x_2\)-axis. Moreover, \(\partial _D I \cup \bigcup _{i=1}^2 \partial _i I\) denotes the subset of boundary points of I that model components of the boundary of \({\mathscr {S}}(t)\). Here endpoints in \(\partial _D I\) correspond to fixed boundary circles of \({\mathscr {S}}(t)\), that lie within a hyperplane parallel to the \(x_1-x_3\)-plane \({{\mathbb {R}}}\times \{0\} \times {{\mathbb {R}}}\). Endpoints in \(\partial _1 I\) correspond to boundary circles of \({\mathscr {S}}(t)\) that can move freely along the boundary of an infinite cylinder that is aligned with the axis of rotation. Endpoints in \(\partial _2 I\) correspond to boundary circles of \({\mathscr {S}}(t)\) that can expand/shrink freely within a hyperplane parallel to the \(x_1-x_3\)-plane \({{\mathbb {R}}}\times \{0\} \times {{\mathbb {R}}}\). See Table 1 for a visualization of the different types of boundary nodes.
On assuming that
we introduce the arclength s of the curve, i.e. \(\partial _s = |\vec {x}_\rho |^{-1}\,\partial _\rho \), and set
where \((\cdot )^\perp \) denotes a clockwise rotation by \(\frac{\pi }{2}\).
On recalling (2.1), we observe that the normal \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) on \({\mathscr {S}}(t)\) is given by
Similarly, the normal velocity \({\mathscr {V}}_{{\mathscr {S}}}\) of \({\mathscr {S}}(t)\) in the direction \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) is given by
For the curvature \(\varkappa \) of \(\varGamma (t)\) it holds that
An important role in this paper is played by the surface area of the surface \({\mathscr {S}}(t)\), which is equal to
Often the surface area, \(A(\vec {x}(t))\), will play the role of the free energy in our paper. But for an open surface \({\mathscr {S}}(t)\), with boundary \(\partial {\mathscr {S}}(t)\), we consider contact energy contributions which are discussed in [24], see also [9, (2.21)]. In the axisymmetric setting the relevant energy is given by
where we recall from (2.2c) that, for \(i=1,2\), either \(\partial _i I = \emptyset \), \(\{0\}\), \(\{1\}\) or \(\{0,1\}\). In the above \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)} \in {{\mathbb {R}}}\), for \(p\in \{0,1\}\), are given constants. Here \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\), for \(p \in \partial _1 I\), denotes the change in contact energy density in the direction of \(-\,\vec {e}_2\), that the two phases separated by the interface \({\mathscr {S}}(t)\) have with the infinite cylinder at the boundary circle of \({\mathscr {S}}(t)\) represented by \(\vec {x}(p,t)\). Similarly, \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\), for \(p \in \partial _2 I\), denotes the change in contact energy density in the direction of \(-\,\vec {e}_1\), that the two phases separated by the interface \({\mathscr {S}}(t)\) have with the hyperplane \({{\mathbb {R}}}\times \{0\}\times {{\mathbb {R}}}\) at the boundary circle of \({\mathscr {S}}(t)\) represented by \(\vec {x}(p,t)\). These changes in contact energy lead to the contact angle conditions
for all \(t \in (0,T]\). In most cases, the contact energies are assumed to be the same, so that \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(0)}={\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)}=0\), which leads to \(90^\circ \) 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 \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)} \not =0 \) for \(p\in \partial _1 I\) or if \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}<0\) for \(p\in \partial _2 I\).
For later use we note that
Moreover, we recall that expressions for the mean curvature and the Gauss curvature of \({\mathscr {S}}(t)\) are given by
respectively; see e.g. [16, (6)]. More precisely, if \(k_m\) and \(k_g\) denote the mean and Gauss curvatures of \({\mathscr {S}}(t)\), then
In the literature, the two terms making up \(\varkappa _{{\mathscr {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
It follows from (2.14) and (2.4) that
A weak formulation of (2.15) 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.17) in the context of a weak formulation of (2.14) will be given in the “Appendix A”.
2.1 Mean curvature flow
In terms of the axisymmetric description of \({\mathscr {S}}(t)\), the evolution law (1.1) can be written as
with, on recalling (2.2b–d),
as well as (2.17) and (2.10a, b).
Let
Then we consider the following weak formulation of (2.18a, b), on recalling (2.7).
\(({\mathscr {A}})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\varkappa (t)\in L^2(I)\) such that
We note that (2.19b) weakly imposes (2.17) and (2.10a, b). We observe that (2.18a) degenerates for \(\vec {x}\cdot \,\vec {e}_1 = 0\), i.e. when \(\rho \in \partial _0 I\). Hence this degeneracy is balanced by the condition (2.16). In fact, on recalling (2.7) it holds that
We remark that the weak formulation \(({\mathscr {A}})\) is close in spirit to the weak formulations introduced in [5, 7] for mean curvature flow. In particular, the tangential component of \(\vec {x}_t\) is not prescribed, which on the discrete level leads to an equidistribution property.
Choosing \(\vec {\eta } = (\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_t \in \underline{V}_\partial \) in (2.19b) and \(\chi = (\vec {x}\cdot \,\vec {e}_1)\,(\vec {x}_t\cdot \,\vec {\nu })\) in (2.19a), we obtain on recalling (2.12), \(\vec {x}_t \in \underline{V}_\partial \), (2.4) and (2.2a) that
An alternative strong formulation of mean curvature flow, in the axisymmetric setting, to (2.18a) is given by
with (2.18b), where we have recalled (2.7). We consider the following weak formulation of (2.22).
\(({\mathscr {B}})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\vec {\varkappa }(t)\in [L^2(I)]^2\) such that
Similarly to (2.19b), we observe that (2.23b) weakly imposes (2.17) and (2.10a, b). We remark that the weak formulation \(({\mathscr {B}})\) in some sense is close in spirit to the weak formulations introduced in [20, 21] for mean curvature flow. In particular, the tangential component of \(\vec {x}_t\) is fixed to be zero, as the right hand side of (2.23a) is normal, recall (2.7).
Choosing \(\vec {\chi } = \vec {\eta } = (\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_t \in \underline{V}_\partial \) in (2.23a, b), we obtain, similarly to (2.21), that
We remark that it does not appear possible to mimic either (2.21) for \(({\mathscr {A}})\) or (2.24) for \(({\mathscr {B}})\) on the discrete level. Hence, in order to develop stable approximations, we investigate alternative formulations based on (2.15). The first formulation corresponds to the strong formulation \((\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_t\cdot \,\vec {\nu } = \vec {x}\cdot \,\vec {e}_1\,\varkappa _{{\mathscr {S}}}\), together with (2.15).
\(({\mathscr {C}})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\varkappa _{{\mathscr {S}}}(t) \in L^2(I)\) such that
The second formulation corresponds to the strong formulation \((\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_t = \vec {x}\cdot \,\vec {e}_1\,\vec {\varkappa }_{{\mathscr {S}}}\), where \(\vec {\varkappa }_{{\mathscr {S}}} = \varkappa _{{\mathscr {S}}}\,\vec {\nu }\), together with (2.15).
\(({\mathscr {D}})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\vec {\varkappa }_{{\mathscr {S}}}(t) \in [L^2(I)]^2\) such that
We note that the variational formulation for \(\vec {\varkappa }_{{\mathscr {S}}}\) in (2.26b) has previously been employed in [26, p. 124].
Choosing \(\vec {\eta } = \vec {x}_t \in \underline{V}_\partial \) in (2.25b) and \(\chi = \varkappa _{{\mathscr {S}}}\) in (2.25a), we obtain for the formulation \(({\mathscr {C}})\), on recalling (2.12), that
Similarly, choosing \(\vec {\eta } = \vec {x}_t \in \underline{V}_\partial \) in (2.26b) and \(\vec {\chi } = \vec {\varkappa }_{{\mathscr {S}}}\) in (2.26a), we obtain for the formulation \(({\mathscr {D}})\) that
We observe that (2.25b) and (2.26b) weakly impose (2.10a, b). But, in contrast to (2.19b) and (2.23b), it is not obvious that they also weakly impose (2.17), due to the presence of the degenerate weight \(\vec {x}\cdot \,\vec {e}_1\). However, we show in the “Appendix A” that in fact they also weakly impose (2.17).
We also note that the formulation \(({\mathscr {C}})\) is loosely related to \(({\mathscr {A}})\), in the sense that the tangential component of \(\vec {x}_t\) is not prescribed. But in contrast to \(({\mathscr {A}})\), discretizations of \(({\mathscr {C}})\) cannot be shown to have an equidistribution property. In a similar way, the formulation \(({\mathscr {D}})\) is loosely related to \(({\mathscr {B}})\), in the sense that the velocity \(\vec {x}_t\) is purely in the normal direction, recall (2.15) and (2.7). Finally, we observe that the variable \(\varkappa _{{\mathscr {S}}}\) can be eliminated from \(({\mathscr {C}})\), by choosing \(\chi = \vec {\nu }\cdot \,\vec {\eta }\) in (2.25a) for \(\vec {\eta } \in \underline{V}_\partial \), and then combining (2.25a) and (2.25b). Similarly, \(\vec {\varkappa }_{{\mathscr {S}}}\) can be eliminated from \(({\mathscr {D}})\) by choosing \(\vec {\chi } = \vec {\eta }\) in (2.26a) for \(\vec {\eta } \in \underline{V}_\partial \), and then combining (2.26a) and (2.26b). We remark that the formulation \(({\mathscr {C}})\), with the variable \(\varkappa _{{\mathscr {S}}}\), as well as the formulation \(({\mathscr {A}})\), are useful with a view towards introducing numerical approximations of higher order flows, such as surface diffusion, see [11].
2.2 Nonlinear mean curvature flow
It is a simple matter to extend the formulations \(({\mathscr {A}})\) and \(({\mathscr {C}})\) to the nonlinear flow (1.2). In principle this can also be achieved for \(({\mathscr {B}})\) and \(({\mathscr {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 \(({\mathscr {A}})\) and \(({\mathscr {C}})\). For the former, replacing the right hand side in (2.19a) with \(\int _I f(\varkappa - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1} )\, \chi \,|\vec {x}_\rho | \;{\mathrm{d}}\rho \) yields a weak formulation for (1.2), which we call \(({\mathscr {A}}^f)\). Similarly, replacing \(\varkappa _{{\mathscr {S}}}\) with \(f(\varkappa _{{\mathscr {S}}})\) in (2.25a) generalizes \(({\mathscr {C}})\) to \(({\mathscr {C}}^f)\) for (1.2).
Similarly to (2.27), and using the same choices of \(\vec {\eta }\) and \(\chi \), it can be shown that solutions to \(({\mathscr {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).
\(({\mathscr {A}}^{f,V})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\varkappa (t)\in L^2(I)\) such that
Choosing \(\chi = 2\,\pi \,\vec {x}\cdot \,\vec {e}_1\) in (2.30a) yields, on recalling (1.5), that
where \({\mathscr {S}}(t) = \partial \varOmega (t)\), and where the sign in (2.31) depends on whether \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) is the outer or inner normal to \(\varOmega (t)\) on \({\mathscr {S}}(t)\), recall (2.5) and (2.6).
\(({\mathscr {C}}^{f,V})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\varkappa _{{\mathscr {S}}}(t) \in L^2(I)\) such that
Choosing \(\chi = 2\,\pi \) in (2.32a) yields (2.31), as before. Moreover, and similarly to (2.29), it can be shown for solutions of \(({\mathscr {C}}^{f,V})\) in the case (1.4c) that
where we have noted the Cauchy–Schwarz inequality. It does not appear possible to extend the stability result (2.33) to the case of more general f.
3 Semidiscrete schemes
Let \([0,1]=\bigcup _{j=1}^J I_j\), \(J\ge 3\), 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={{\mathbb {R}}} \slash {{\mathbb {Z}}}\) we identify \(0=q_0 = q_J=1\).
The necessary finite element spaces are given by
\(\underline{V}^h= [V^h]^2\), \(\underline{V}^h_{\partial _0}= \underline{V}^h\cap \underline{V}_{\partial _0}\) and \(\underline{V}^h_\partial = \underline{V}^h\cap \underline{V}_\partial \). Let \(\{\chi _j\}_{j=j_0}^J\) denote the standard basis of \(V^h\), where \(j_0 = 0\) if \(I = (0,1)\) and \(j_0 = 1\) if \(I={{\mathbb {R}}} \slash {{\mathbb {Z}}}\). For later use, we let \(\pi ^h{:}\,C({\overline{I}})\rightarrow V^h\) be the standard interpolation operator at the nodes \(\{q_j\}_{j=0}^J\). We extend \(\pi ^h\) to \(\vec \pi ^h\) for vector valued functions. Let \((\cdot ,\cdot )\) 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=1}^J\), via
where we define \(f(q_j^\pm )=\lim \nolimits _{\delta \searrow 0}f(q_j\pm \delta )\). The definition (3.2) naturally extends to vector valued functions. It is easily shown that
Let \((\vec {X}^h(t))_{t\in [0,T]}\), with \(\vec {X}^h(t)\in \underline{V}^h_{\partial _0}\), be an approximation to \((\vec {x}(t))_{t\in [0,T]}\) and define \(\varGamma ^h(t) = \vec {X}^h(t)({\overline{I}})\). Throughout this section we assume that
Assuming that \(|\vec {X}^h_\rho | > 0\) almost everywhere on I, and similarly to (2.4), we set
We note that
For later use, we let \(\vec {\omega }^h \in {\underline{V}}^h\) be the mass-lumped \(L^2\)-projection of \(\vec {\nu }^h\) onto \({\underline{V}}^h\), i.e.
Recall from (2.8) and (2.9) that
We have, similarly to (2.12), that
3.1 Mean curvature flow
In view of the degeneracy on the right hand side of (2.18a), and on recalling (2.20) and (3.6), we introduce, given a \(\kappa ^h(t) \in V^h\), the function \({\mathfrak {K}}^h(\kappa ^h(t),t) \in V^h\) such that
Our semidiscrete finite element approximation of \(({\mathscr {A}})\), (2.19a, b), is given as follows.
\(({\mathscr {A}}_h)^h\) Let \(\vec {X}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {X}^h(t) \in \underline{V}^h\), with \(\vec {X}^h_t(t) \in \underline{V}^h_\partial \), and \(\kappa ^h(t) \in V^h\) such that
Remark 3.1
Let \(\vec {h}_j(t) = \vec {X}^h(q_j,t) - \vec {X}^h(q_{j-1},t)\) for \(j = 1,\ldots , J\), and set \(\vec {h}_0 = \vec {h}_J\) if \(\partial I = \emptyset \). Then, if \((\vec {X}^h(t), \kappa ^h(t)) \in \underline{V}^h\times V^h\) satisfies (3.10b), it holds that
The equidistribution property (3.11) can be shown by choosing \(\vec {\eta }=\chi _{j-1}\,[\vec {\omega }^h(q_{j-1},t)]^\perp \in \underline{V}^h_\partial \) in (3.10b), recall (3.6). See also [6, Remark 2.4] and [5, Remark 2.5] for more details.
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 \(({\mathscr {A}}_h)^h\) with mass lumping. Of course, in the case \(\partial _0 I= \emptyset \), an alternative scheme to (3.10a, b) is
together with (3.10b). Note that if \(\partial _0 I= \emptyset \) then (3.10a) collapses to (3.13) with \(\vec {\nu }^h\) replaced by \(\vec {\omega }^h\). Unfortunately, neither choice appears to lead to a stability proof.
In an attempt to prove stability, we choose \(\vec {\eta } = \vec {\pi }^h[(\vec {X}^h\cdot \,\vec {e}_1)\,\vec {X}^h_t]\) in (3.10b). Then it follows from (3.8), \(\vec {X}^h_t \in \underline{V}^h_\partial \), (3.5) and (3.6) that
Moreover, considering for simplicity the case \(\partial _0 I= \emptyset \), and choosing \(\chi = -\pi ^h[2\,\pi \,(\vec {X}^h\cdot \,\vec {e}_1)\,(\vec {X}^h_t\cdot \,\vec {\omega }^h)]\) in (3.13) yields, on noting (3.6), that
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 \((\vec {X}^h\cdot \,\vec {e}_1)\,(\vec {X}^h_t\cdot \,\vec {\nu }^h)\) is discontinuous, and so \(\pi ^h[(\vec {X}^h\cdot \,\vec {e}_1)\,(\vec {X}^h_t\cdot \,\vec {\nu }^h)]\) is not well-defined, and cannot be chosen as a test function in (3.10a) or (3.13).
However, the fully discrete variant of \(({\mathscr {A}}_h)^h\), (3.10a, b), in our numerical simulations did not exhibit any stability issues.
A semidiscrete approximation of \(({\mathscr {B}})\), (2.23a, b), is given as follows.
\(({\mathscr {B}}_h)^h\) Let \(\vec {X}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {X}^h(t) \in \underline{V}^h\), with \(\vec {X}^h_t(t) \in \underline{V}^h_\partial \), and \(\vec {\kappa }^h(t) \in \underline{V}^h\), such that
where \(\vec {{\mathfrak {K}}}^{h}(\vec {\kappa }^h) \in \underline{V}^h\) is such that
The rescaling factor \(|\vec {\omega }^h(q_j,t)|^2\) in (3.17) normalizes the discrete vertex normals \(\vec {\omega }^h(q_j,t)\), recall (3.12), which is the most natural approach. Similarly to \(({\mathscr {A}}_h)^h\), it does not appear possible to prove a stability result for \(({\mathscr {B}}_h)^h\).
However, it turns out that approximations of the formulations \(({\mathscr {C}})\) and \(({\mathscr {D}})\) can be shown to be stable. In particular, our semidiscrete approximations of \(({\mathscr {C}})\), (2.25a, b), and \(({\mathscr {D}})\), (2.26a, b), are given as follows, where we first define
\(({\mathscr {C}}_h)^{(h)}\) Let \(\vec {X}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {X}^h(t) \in \underline{V}^h\), with \(\vec {X}^h_t(t) \in \underline{V}^h_\partial \), and \(\kappa _{{\mathscr {S}}}^h(t) \in W^h_{(\partial _0)}\) such that
Here and throughout we use the notation \(\cdot ^{(h)}\) to denote an expression with or without the superscript h, and similarly for the subscripts \(\cdot _{(\partial _0)}\). I.e. the scheme \(({\mathscr {C}}_h)^h\) employs mass lumping, recall (3.2), and seeks \(\kappa _S(t) \in W^h_{\partial _0}\), while the scheme \(({\mathscr {C}}_h)\) employs true integration throughout and seeks \(\kappa _S(t) \in W^h = V^h\).
\(({\mathscr {D}}_h)^{(h)}\) Let \(\vec {X}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {X}^h(t) \in \underline{V}^h\), with \(\vec {X}^h_t(t) \in \underline{V}^h_\partial \), and \(\vec {\kappa }_{{\mathscr {S}}}^h(t) \in \underline{W}_{(\partial _0)}^h\) such that
We observe that \(({\mathscr {C}}_h)^{h}\) and \(({\mathscr {D}}_h)^{h}\) do not depend on the values of \(\kappa _{{\mathscr {S}}}^h\) and \(\vec {\kappa }_{{\mathscr {S}}}^h\), respectively, on \(\partial _0 I\). Hence we fix these values to be zero by requiring that \(\kappa _{{\mathscr {S}}}^h \in W^h_{\partial _0}\) and \(\vec {\kappa }_{{\mathscr {S}}}^h \in {\underline{W}}^h_{\partial _0}\), and by using a reduced set of test functions in (3.18a) and (3.19a). As a consequence, it seems at first that \(\vec {X}^h_t\) is not defined on \(\partial _0 I\). However, \(\vec {X}^h\) on \(\partial _0 I\) is determined through (3.18b) and (3.19b), respectively.
We have on choosing \(\chi = \kappa _{{\mathscr {S}}}^h\) in (3.18a), \(\vec {\chi } = \vec {\kappa }_{{\mathscr {S}}}^h\) in (3.19a) and \(\vec {\eta } = \vec {X}^h_t\) in (3.18b) and (3.19b), on recalling (3.8), that
respectively. This shows that both methods are stable, where we recall (3.4). Similarly to (2.28) and (2.27), we observe that (3.20) implies for \(({\mathscr {D}}_h)^{(h)}\) and \(({\mathscr {C}}_h)^h\) 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 \(\vec {\kappa }^h_{{\mathscr {S}}}\) from the schemes \(({\mathscr {D}}_h)^{(h)}\), which yields (3.19b) with \(\vec {\kappa }^h_{{\mathscr {S}}}\) replaced by \(\vec {X}^h_t\). Similarly, \(\kappa ^h_{{\mathscr {S}}}\) can be removed from the scheme \(({\mathscr {C}}_h)^{h}\) to yield (3.18b) with \(\kappa ^h_{{\mathscr {S}}}\,\vec {\nu }^h\) replaced by \((\vec {X}^h_t\cdot \,\vec {\omega }^h)\,\vec {\omega }^h\), on recalling (3.6). For the scheme \(({\mathscr {C}}_h)\) this elimination procedure is not possible.
For the reader’s convenience, Table 2 summarises the main properties of all the schemes introduced in this section.
3.2 Nonlinear mean curvature flow
Replacing \(\kappa ^h - {\mathfrak {K}}^h(\kappa ^h)\) with \(f(\kappa ^h - {\mathfrak {K}}^h(\kappa ^h))\) in (3.10a) yields the scheme \(({\mathscr {A}}_h^f)^h\). Similarly, the scheme \(({\mathscr {A}}_h^{f,V})^h\) is given by (3.10a, b) with the right hand side in (3.10a) replaced by
These two schemes inherit the equidistribution property, recall (3.11). Replacing \(\kappa ^h_{{\mathscr {S}}}\) with \(\pi ^h[f(\kappa ^h_{{\mathscr {S}}})]\) in (3.18a) yields the schemes \(({\mathscr {C}}_h^f)^{(h)}\) and similarly we can define \(({\mathscr {C}}_h^{f,V})^{(h)}\) by replacing the right hand side in (3.18a) by
Similarly to (3.20), and using the same choices of \(\vec {\eta }\) and \(\chi \), it can be shown that solutions to the scheme \(({\mathscr {C}}_h^f)^{(h)}\) satisfy \(-\frac{1}{2\,\pi }\,\frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {X}^h(t)) = \left( (\vec {X}^h\cdot \,\vec {e}_1)\,f(\kappa _{{\mathscr {S}}}^h), \kappa _{{\mathscr {S}}}^h\, |\vec {X}^h_\rho |\right) ^{(h)}\), which yields a stability bound for \(({\mathscr {C}}_h^f)^{h}\) if f is monotonically increasing with \(f(0) = 0\). Of course, (3.22) is a discrete analogue of (2.29). Moreover, solutions to \(({\mathscr {C}}_h^{f,V})^{(h)}\), in the case (1.4c), satisfy
similarly to (2.33), where here we have also used a Cauchy–Schwarz inequality for the mass lumped inner product (3.2). Finally, solutions to the scheme \(({\mathscr {C}}_h^{f,V})\) conserve the volume of the domain \(\varOmega ^h(t) \subset {{\mathbb {R}}}^3\) that is enclosed by the three-dimensional axisymmetric surface \({\mathscr {S}}^h(t)\) that is generated by the curve \(\varGamma ^h(t)\). To see this, choose \(\chi = 2\,\pi \) in (3.18a), with the modified right hand side (3.22), to obtain
recall (1.5). Here \({\mathscr {V}}^h_{{\mathscr {S}}^h}(t)\) denotes the normal velocity of \({\mathscr {S}}^h(t)\) in the direction of \(\vec {\nu }^h_{{\mathscr {S}}^h}(t)\), the outer normal to \(\varOmega ^h(t)\) on \({\mathscr {S}}^h(t)\), where \(\vec {\nu }^h_{{\mathscr {S}}^h}(t)\) is induced by \(\vec {\nu }^h\) through a discrete analogue of (2.5). Using the same testing procedure for the scheme \(({\mathscr {C}}_h^{f,V})^h\) yields that
and so the enclosed volume is only approximately preserved, compare with (3.23). Finally, choosing \(\chi = \vec {X}^h\cdot \,\vec {e}_1\) in \(({\mathscr {A}}_h^{f,V})^h\), recall (3.21), also yields (3.24), and so an approximate volume preservation property.
4 Fully discrete schemes
In this section we present and analyse fully discrete variants of the approximations discussed in Sect. 3 that turn out to be most useful in practice. Other fully discrete schemes, which we will use for comparison purposes in the numerical results section, are summarized in “Appendix C”.
Let \(0= t_0< t_1< \cdots< t_{M-1} < t_M = T\) be a partitioning of [0, T] into possibly variable time steps \(\varDelta t_m = t_{m+1} - t_{m}\), \(m=0\rightarrow M-1\). We set \(\varDelta t= \max _{m=0\rightarrow M-1}\varDelta t_m\). For a given \(\vec {X}^m\in \underline{V}^h_{\partial _0}\), assuming that \(|\vec {X}^m_\rho | > 0\) almost everywhere on I, we set \(\vec {\nu }^m = - \frac{[\vec {X}^m_\rho ]^\perp }{|\vec {X}^m_\rho |}\). Let \(\vec {\omega }^m \in \underline{V}^h\) be the natural fully discrete analogue of \(\vec {\omega }^h \in \underline{V}^h\), recall (3.6), i.e.
4.1 Mean curvature flow
Similarly to (3.9), and given a \(\kappa ^{m+1} \in V^h\), we introduce \({\mathfrak {K}}^{m}(\kappa ^{m+1}) \in V^h\) such that
Then our fully discrete analogue of \(({\mathscr {A}}_h)^h\), (3.10a, b), is given as follows.
\(({\mathscr {A}}_m)^h\) Let \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \((\delta \vec {X}^{m+1}, \kappa ^{m+1}) \in \underline{V}^h_\partial \times V^h\), where \(\vec {X}^{m+1} = \vec {X}^m + \delta \vec {X}^{m+1}\), such that
We make the following assumptions.
Note that the assumption \(({\mathfrak {B}})^h\), on recalling (3.6), is equivalent to assuming that \(\dim \text {span}\{\vec {\omega }^m(q_j)\}_{j=0}^{J}= 2\). Clearly, it can only be violated if all the \(J+1\) vectors lie on a single line. It was shown in [6, Remark 2.2] that the assumption \(({\mathfrak {B}})^h\) always holds for closed curves without self-intersections.
Lemma 4.1
Let \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfy the assumptions \(({\mathfrak {A}})\) and \(({\mathfrak {B}})^h\). Then there exists a unique solution \((\delta \vec {X}^{m+1}, \kappa ^{m+1}) \in \underline{V}^h_\partial \times V^h\) to \(({\mathscr {A}}_m)^h\).
Proof
We note that since \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfies the assumption \(({\mathfrak {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: find \((\delta \vec {X},\kappa ) \in \underline{V}^h_\partial \times V^h\) such that
where we recall from (4.2) that \(\lambda \in V^h\) with
Choosing \(\chi =\kappa \in V^h\) in (4.4a) and \(\vec {\eta }= \delta \vec {X} \in \underline{V}^h_\partial \) in (4.4b) yields that
It follows from (4.6) that \(\kappa = 0\) and that \(\delta \vec {X} \equiv \vec {X}^c\in {{\mathbb {R}}}^2\); and hence that
It follows from (4.7) and assumption \(({\mathfrak {B}})^h\) that \(\vec {X}^c=\vec {0}\). Hence we have shown that (4.3a, b) has a unique solution \((\delta \vec {X}^{m+1},\kappa ^{m+1}) \in \underline{V}^h_\partial \times V^h\).
A fully discrete analogue of the scheme \(({\mathscr {C}}_h)^{(h)}\), (3.18a, b), is given as follows. Here we introduce the notation \([ r ]_\pm = \pm \max \{ \pm r, 0 \}\) for \(r \in {{\mathbb {R}}}\).
\(({\mathscr {C}}_{m,\star })^{(h)}\) Let \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \((\delta \vec {X}^{m+1}, \kappa _{{\mathscr {S}}}^{m+1}) \in \underline{V}^h_\partial \times W^h_{(\partial _0)}\), where \(\vec {X}^{m+1} = \vec {X}^m + \delta \vec {X}^{m+1}\), such that
We note that in practice we solve the nonlinear scheme \(({\mathscr {C}}_{m,\star })^{(h)}\) with a Newton iteration, which in all our experiments always converged with at most three iterations.
Remark 4.1
Similarly to the semidiscrete variants, we observe that it is possible to eliminate the discrete curvature, \(\kappa _{{\mathscr {S}}}^{m+1}\) from \(({\mathscr {C}}_{m,\star })^h\). On recalling (3.6) and on choosing \(\chi = \pi ^h[\vec {\eta }\cdot \,\vec {\omega }^m] \in W^h_{(\partial _0)}\) in (4.8a) for \(\vec {\eta } \in \underline{V}^h_\partial \), the scheme reduces to: Find \(\delta \vec {X}^{m+1}\in \underline{V}^h_\partial \) such that, for all \(\vec {\eta } \in \underline{V}^h_\partial \),
For the scheme \(({\mathscr {C}}_{m,\star })\) this elimination procedure is not possible. A related variant to (4.9) is given by: find \(\delta \vec {X}^{m+1}\in \underline{V}^h_\partial \) such that, for all \(\vec {\eta } \in \underline{V}^h_\partial \),
Theorem 4.1
Let \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfy the assumption \((\mathfrak A)\). Let \((\vec {X}^{m+1},\kappa _{{\mathscr {S}}}^{m+1})\) be a solution to \(({\mathscr {C}}_{m,\star })^{(h)}\). Then it holds that
where we recall the definition (3.7).
Proof
Choosing \(\chi = \varDelta t_m\,\kappa _{{\mathscr {S}}}^{m+1}\) in (4.8a) and \(\vec {\eta } = \vec {X}^{m+1} - \vec {X}^m \in \underline{V}^h_\partial \) in (4.8b) yields, on noting that \(\vec {X}^{m}(p)\cdot \,\vec {e}_1 = \vec {X}^{m+1}(p)\cdot \,\vec {e}_1\) for \(p \in \partial _1 I\), that
where we have used the two inequalities \(\vec {a}\cdot \,(\vec {a} - \vec {b}) \ge |\vec {b}|\,(|\vec {a}| - |\vec {b}|)\) for \(\vec {a}\), \(\vec {b} \in {{\mathbb {R}}}^2\), and \(2\,\gamma \,(\gamma - \alpha ) \ge \gamma ^2 - \alpha ^2\) for \(\alpha ,\gamma \in {{\mathbb {R}}}\). This proves the desired result.
We note that the scheme (4.10) can also be shown to be unconditionally stable, i.e. a solution to (4.10) satisfies
4.2 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 \(({\mathscr {A}}^f_{m})^{h}\) of the scheme \(({\mathscr {A}}_{m})^{h}\) to (1.2) is given by (4.3a, b) with (4.3a) replaced by
Similarly, the natural adaptation \(({\mathscr {C}}^f_{m,\star })^{(h)}\) of the scheme \(({\mathscr {C}}_{m,\star })^{(h)}\) to (1.2) is given by (4.8a, b) with (4.8a) replaced by
Similarly to (4.11), with the same choices of \(\chi \) and \(\vec {\eta }\), it is then possible to prove that solutions to \(({\mathscr {C}}^f_{m,\star })^{h}\) satisfy \(E(\vec {X}^{m+1}) + 2\,\pi \,\varDelta t_m\, \left( \vec {X}^m\cdot \,\vec {e}_1\, f( \kappa _{{\mathscr {S}}}^{m+1}), \kappa _{{\mathscr {S}}}^{m+1} \,|\vec {X}^m_\rho | \right) ^{h} \le E(\vec {X}^m)\), which provides a stability bound if f is monotonically increasing with \(f(0)=0\).
Finally, replacing (4.13) with
and replacing (4.14) with
gives the fully discrete approximations \(({\mathscr {A}}^{f,V}_{m})^{h}\) and \(({\mathscr {C}}^{f,V}_{m,\star })^{(h)}\), respectively, of (1.3). For the case (1.4c) it is possible to prove a stability bound for \(({\mathscr {C}}^{f,V}_{m,\star })^{(h)}\). In particular, solutions to \(({\mathscr {C}}^{f,V}_{m,\star })^{(h)}\) satisfy, similarly to (4.12), that
where we have noted the Cauchy–Schwarz inequality. For the fully discrete approximations \(({\mathscr {A}}^{f,V}_{m})^{h}\) and \(({\mathscr {C}}^{f,V}_{m,\star })^{(h)}\) it is not possible to prove a volume conservation property. However, in practice all three schemes preserve the enclosed volume well, with the relative volume loss decreasing as the discretization parameters become smaller.
Finally, we note that for the schemes \(({\mathscr {A}}_m^f)^h\) and \(({\mathscr {A}}_m^{f,V})^h\), depending on the choice of f, existence and uniqueness results can be shown, see “Appendix B”.
Remark 4.2
In order to be able to compute evolutions for the general flow (1.6), we propose the scheme \(({\mathscr {A}}^F_{m})^{h}\), which can be obtained from the scheme \(({\mathscr {A}}_{m})^{h}\) (4.3a, 4.3b), by replacing (4.3a) with
This is linear scheme for which existence of a unique solution, provided that the assumption \(({\mathfrak {B}})^h\) holds, can easily be shown. Moreover, solutions to the semidiscrete variant of \(({\mathscr {A}}^F_{m})^{h}\) satisfy the equidistribution property (3.11).
5 Numerical results
As the fully discrete energy, we consider \(E(\vec {X}^m)\), recall (3.7). We always employ uniform time steps, \(\varDelta t_m = \varDelta t\), \(m=0,\ldots ,M-1\).
We also consider the ratio
between the longest and shortest element of \(\varGamma ^m\), and are often interested in the evolution of this ratio over time.
On recalling (2.7), and given \(\varGamma ^0 = \vec {X}^0({\overline{I}})\), for the scheme \(({\mathscr {A}}_m^{f,V})^h\) we define \(\kappa ^0 \in V^h\) via
recall (4.1), where \(\vec {\kappa }^0\in \underline{V}^h\) is such that
5.1 Numerical results for mean curvature flow
5.1.1 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 \(\partial _0 I = \partial I = \{0,1\}\) and we choose \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\) with
recall (3.1). We compute the error
over the time interval [0, T] between the true solution (5.2) and the discrete solutions for the schemes \(({\mathscr {A}}_m)^h\), \(({\mathscr {B}}_m)^h\), \(({\mathscr {C}}_m)^{(h)}\), \(({\mathscr {D}}_m)^{(h)}\), \(({\mathscr {C}}_{m,\star })^{(h)}\) and \(({\mathscr {D}}_{m,\star })^{(h)}\). Here we used the time step size \(\varDelta t=0.1\,h^2_{\varGamma ^0}\), where \(h_{\varGamma ^0}\) is the maximal edge length of \(\varGamma ^0\). The computed errors are reported in Tables 3, 4 and 5. Comparing the reported numbers with the values in 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 3 that the schemes \(({\mathscr {A}}_m)^h\) and \(({\mathscr {B}}_m)^h\) appear to converge with the optimal convergence rate of \({\mathscr {O}}(h^2_{\varGamma ^0})\). Similarly, Tables 4 and 5 suggest that the schemes \(({\mathscr {C}}_{m,\star })^{(h)}\), \(({\mathscr {D}}_m)^{(h)}\) and \(({\mathscr {D}}_{m,\star })^{(h)}\) converge with an order slightly less than quadratic. We note that the linear schemes \(({\mathscr {C}}_m)^{(h)}\) lead to solutions \(\vec {X}^m\) with \(\min _{\rho \in {\overline{I}}} \vec {X}^m(\rho )\cdot \,\vec {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 \(({\mathscr {D}}_m)^{(h)}\) behave well, on the other hand. But as they are very close to the nonlinear schemes \(({\mathscr {D}}_{m,\star })^{(h)}\), from now on we concentrate on the schemes \(({\mathscr {D}}_{m,\star })^{(h)}\), \(({\mathscr {C}}_{m,\star })^{(h)}\), \(({\mathscr {A}}_{m})^{h}\) and \(({\mathscr {B}}_{m})^{h}\).
5.1.2 Torus
We repeat the two torus experiments in Figures 5 and 6 in [7]. To this end, we let \(\partial I = \emptyset \). For an initial torus with radii \(R=1\), \(r=0.7\), we obtain a surface that closes up towards a genus-0 surface, as in [7, Fig. 5]. See Fig. 2 for the simulation results for the scheme \(({\mathscr {A}}_m)^h\), for the discretization parameters \(J=256\) and \(\varDelta 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 Fig. 3 for the evolution for the scheme \(({\mathscr {A}}_m)^h\), again for \(J=256\) and \(\varDelta t= 10^{-4}\). On repeating the numerical experiment for the \(({\mathscr {B}}_m)^h\), we observe strong oscillations, as shown on the right of Fig. 3. These oscillations become smaller in magnitude as \(\varDelta t\) is decreased. The remaining schemes can integrate the evolution shown in Fig. 3 in a stable way, and their numerical results are very close to the ones displayed in Fig. 3 for the scheme \(({\mathscr {A}}_m)^h\). However, the schemes differ in the exhibited tangential motions, which leads to diverse evolutions of the ratio \({{\mathfrak {r}}}^m\), see Fig. 4. The best distribution of mesh points is shown by the schemes \(({\mathscr {A}}_m)^h\) and \(({\mathscr {C}}_{m,\star })^h\), followed by \(({\mathscr {C}}_{m,\star })\). The most nonuniform distribution of mesh points can be observed for the two schemes \(({\mathscr {D}}_{m,\star })^{(h)}\).
5.1.3 Cylinder
For the scheme \(({\mathscr {A}}_m)^h\) we repeat the singular evolution from [22, Fig. 1]. To this end, we set \(\partial _D I = \partial I = \{0,1\}\). In particular, starting with a cylinder, mean curvature flow leads to a pinch-off. We show the results for the scheme \(({\mathscr {A}}_m)^h\), with the discretization parameters \(J=128\) and \(\varDelta t= 10^{-4}\), in Fig. 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 \(\partial _2 I = \partial I = \{0,1\}\) and use the discretization parameters \(J=128\) and \(\varDelta t= 10^{-3}\). Letting \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(0)} = {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)} = -\frac{1}{2}\) and starting with a cylinder, the evolution yields a growing catenoid-like surface, see Fig. 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 [37, p. 15] and [28],
where \(z_0 \in {{\mathbb {R}}}\) specifies the position of the travelling wave solution at time \(t=0\). In fact, plotting \(\vec {g} (\rho ,t) - (z_0 + \tfrac{\pi }{3}\,t)\,\vec {e}_1\) at time \(t=4\) versus \(\vec {X}^m - (\min _{\rho \in {\overline{I}}} \vec {X}^m\cdot \,\vec {e}_1)\,\vec {e}_1\) for our final solution in Fig. 6, yields perfect agreement between the two graphs. We conjecture that the speed of the travelling wave type solution will approach \(\frac{\pi }{3}\), the speed of (5.5). To test this conjecture, we continue the evolution until \(t=100\) and plot the evolution of \(\vec {X}^m(0)\cdot \,\vec {e}_1\) over time, comparing the graph with a suitably chosen line with slope \(\frac{\pi }{3}\), see Fig. 7. As we can see, the speed of the curve does indeed approach \(\frac{\pi }{3}\).
If we let \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(0)} = {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)} = \frac{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 Fig. 8 for the evolution. To allow for an accurate resolution of the detaching, we employ the smaller time step size \(\varDelta t=10^{-6}\) for this simulation.
5.1.4 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 \(\partial _0 I = \{0\}\) and \(\partial _1 I = \{1\}\), and use the discretization parameters \(J=128\) and \(\varDelta t= 10^{-3}\). Letting \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)} = -\frac{1}{2}\) and starting with a disk, the evolution seems to converge to a translating surface patch, see Fig. 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.
5.2 Numerical results for conserved mean curvature flow
5.2.1 Sphere
Clearly, a sphere is a stationary solution for conserved mean curvature flow, (1.3) with (1.4c). Hence, setting \(\partial _0 I = \partial I = \{0,1\}\) and choosing as initial data the nonuniform approximation of a semicircle (5.3) with \(J=64\), we now investigate the different tangential motions exhibited by our proposed schemes. The initial data \(\vec {X}^0\) has a ratio \({{\mathfrak {r}}}^0=1.22\), recall (5.1). We set \(\varDelta t= 10^{-4}\) and integrate the evolution until time \(T=1\). For the three schemes \(({\mathscr {A}}_m)^h\), \(({\mathscr {C}}_{m,\star })^h\) and \(({\mathscr {C}}_{m,\star })\) the element ratios \({{\mathfrak {r}}}^m\) at time \(T=1\) are 1.01, 73.13, 2.94, and the enclosed volume is preserved almost exactly by all the schemes. We show the final distributions of vertices, and plots of \({{\mathfrak {r}}}^m\) over time in Fig. 10. An insight that we gain from this set of experiments is that the tangential motion displayed by the scheme \(({\mathscr {C}}_{m,\star })^h\) can lead to very nonuniform meshes. Hence, for the remainder of this paper, we will only present numerical results for the two schemes \(({\mathscr {A}}_m)^h\) and \(({\mathscr {C}}_{m,\star })\) and their nonlinear variants. Note that the former is a linear fully discrete approximation of \(({\mathscr {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.1. As the results for \(({\mathscr {A}}_m)^h\) and \(({\mathscr {C}}_{m,\star })\) are often indistinguishable, we only visualize the numerical results for the former from now on.
5.2.2 Genus 0 surface
An experiment for a cigar shape can be seen in Fig. 11. Here we have once again that \(\partial _0 I = \partial I = \{0,1\}\). The discretization parameters are \(J=128\) and \(\varDelta t= 10^{-4}\). The relative volume loss for this experiment for \(({\mathscr {A}}_m^{f,V})^h\) is \(0.09\%\), while for \(({\mathscr {C}}_{m,\star }^{f,V})\) it is \(-\,0.01\%\). The same experiment for the scheme \(({\mathscr {C}}_{m,\star }^{f,V})^h\) yields a very nonuniform mesh, with the final ratio \({{\mathfrak {r}}}^m > 430\).
An experiment for a disc shape is shown in Fig. 12. The discretization parameters are \(J=128\) and \(\varDelta t= 10^{-4}\). The relative volume loss for this experiment for \(({\mathscr {A}}_m^{f,V})^h\) is \(-\,0.02\%\), while for \(({\mathscr {C}}_{m,\star }^{f,V})\) it is \(-\,0.01\%\). Once again, the scheme \(({\mathscr {C}}_{m,\star }^{f,V})^h\) yields a very nonuniform mesh for this simulation, with the final ratio \({{\mathfrak {r}}}^m > 145\).
5.2.3 Genus 1 surface
We repeat the simulation in Fig. 3 for conserved mean curvature flow, i.e. (1.3) with (1.4c), using the scheme \(({\mathscr {A}}_m^{f,V})^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 Fig. 13. As for the original experiment, we use the discretization parameters \(J=256\) and \(\varDelta t= 10^{-4}\). The relative enclose volume loss for this experiment is \(-\,0.00\%\). The evolutions for the schemes \(({\mathscr {C}}_{m,\star }^{f,V})^{h}\) and \(({\mathscr {C}}_{m,\star }^{f,V})\) are nearly identical to what is shown in Fig. 13, with a relative volume loss of \(0.01\%\) in both cases.
Finally, we present an example for conserved mean curvature flow, (1.3) with (1.4c), for the scheme \(({\mathscr {A}}_m^{f,V})^h\) with the initial data \(\vec {X}^0\) parameterizing a closed spiral, so that the approximated surface has genus 1. As can be seen from Fig. 14, the spiral slowly untangles, until the surface becomes a torus. For this experiment we use the discretization parameters \(J=1024\) and \(\varDelta t= 10^{-6}\). The relative enclosed volume loss for this experiment is \(0.01\%\). The evolutions for the schemes \(({\mathscr {C}}_{m,\star }^{f,V})^{h}\) and \(({\mathscr {C}}_{m,\star }^{f,V})\) are nearly identical to what is shown in Fig. 14, with a relative volume loss of \(0.01\%\) in both cases.
5.3 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 \(\beta = \tfrac{1}{2}\), similarly to Table 2 in [7]. Here we start with the nonuniform partitioning (5.3) of a semicircle of radius \(r(0)=r_0=1\) and compute the flow until time \(T = \tfrac{1}{2}\,{\overline{T}}\), where \({\overline{T}} = \tfrac{2}{3}\,2^{-\frac{1}{2}}\) denotes that extinction time of the shrinking sphere. We compute the error \(\Vert \varGamma - \varGamma ^h\Vert _{L^\infty }\), recall (5.4), over the time interval [0, T] between the true solution (5.6) and the discrete solutions for the schemes \(({\mathscr {A}}_m^f)^h\) and \(({\mathscr {C}}_{m,\star }^f)\). Here we used the time step size \(\varDelta t=0.1\,h^2_{\varGamma ^0}\), where \(h_{\varGamma ^0}\) is the maximal edge length of \(\varGamma ^0\). The computed errors are reported in Table 6.
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 7. We recall that these numbers can be compared to the fully 3d results in Table 3 in [7]. It is clear from Tables 6 and 7 that the solutions to the scheme \(({\mathscr {A}}_m^f)^h\) appear to converge with the optimal convergence rate of \({\mathscr {O}}(h^2_{\varGamma ^0})\). For the scheme \(({\mathscr {C}}_{m,\star }^f)\), on the other hand, the solutions appear to converge with an order less than quadratic, and closer to \(\frac{3}{2}\). We believe that these lower convergence rates are caused by the nonuniform meshes induced by the scheme \(({\mathscr {C}}_{m,\star }^f)\), recall Fig. 10, and also by the degeneracy of the coefficients \(\vec {x}\cdot \,\vec {e}_1\) in \(({\mathscr {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 \(\varkappa _{{\mathscr {S}}}(\cdot ,0)<0\), the classical inverse mean curvature develops a singularity in finite time, see also [27, 43]. For the axisymmetric setting we use \(I = {{\mathbb {R}}}/{\mathbb {Z}}\), so that \(\partial I = \emptyset \). As the discretization parameters for the scheme \(({\mathscr {A}}_m^f)^h\) we use \(J=256\) and \(\varDelta t= 10^{-4}\). See Fig. 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.
5.4 Numerical results for Gauss curvature flow
An experiment for Gauss curvature flow, (1.7), for the same initial data as in Fig. 11, can be seen in Fig. 16. Here we have once again that \(\partial _0 I = \partial I = \{0,1\}\). The discretization parameters for the scheme \(({\mathscr {A}}^F_{m})^{h}\) from Remark 4.2 are \(J=128\) and \(\varDelta t= 10^{-5}\). As a comparison, we also show the evolution for standard mean curvature flow, computed with the scheme \(({\mathscr {A}}_m)^{h}\), in Fig. 16. It was suggested by Firey [25], that surfaces of stones, which are pounded by waves and other stones, move according to Gauss curvature flow. It is more likely that parts of the surface, where both principal curvature directions are highly curved, will be hit by waves and other stones. He hence proposed the Gauss curvature flow as the governing equation for the evolution of the stone’s surface. In Fig. 16 it is clearly seen that the upper and lower part, which have two highly curved principal curvature directions, move faster within Gauss curvature flow when compared to mean curvature flow. The parts closer to the origin have a nearly flat principal curvature direction. Hence they move far slower under Gauss curvature flow than under mean curvature flow, as can be clearly seen in Fig. 16.
The Gauss curvature flow is not well-defined for general hypersurfaces as for non-convex hypersurfaces the resulting equation is not parabolic, see [3, 34]. In the axisymmetric situation the degrees of freedom are reduced and the resulting equation is
which is parabolic as long as \({\vec {\nu }\cdot \,\vec {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 Fig. 17, where we used the discretization parameters \(J=128\) and \(\varDelta t= 10^{-5}\). For a mathematical analysis for Gauss curvature flow in the axisymmetric case we refer to [34].
6 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 \(({\mathscr {A}}_m)^h\) and \(({\mathscr {C}}_{m,\star })\). Here we have dismissed the scheme \(({\mathscr {B}}_m)^h\), as it can show oscillations in practice, recall Fig. 3, as well as the scheme \(({\mathscr {C}}_{m,\star })^h\), as it can display very nonuniform meshes in simulations where the discrete curves are attached to the \(x_2\)-axis, recall Fig. 10 for its variant \(({\mathscr {C}}_{m,\star }^{f,V})^h\). We also do not consider the schemes \(({\mathscr {D}}_{m,\star })^{(h)}\), as they have no advantage over \(({\mathscr {C}}_{m,\star })\) and as they can also exhibit very nonuniform meshes. Of the two schemes we consider, the scheme \(({\mathscr {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 \(({\mathscr {A}}_m)^h\), in practice the discrete energy is always monotonically decreasing. The scheme \(({\mathscr {C}}_{m,\star })\), 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 \(({\mathscr {C}}_{m,\star })\) may be worse than for \(({\mathscr {A}}_m)^h\), but coalescence of vertices is not observed in practice. Similar statements hold for the nonlinear variants \(({\mathscr {A}}_m^f)^h\), \(({\mathscr {A}}_m^{f,V})^h\), \(({\mathscr {C}}_{m,\star }^f)\) and \(({\mathscr {C}}_{m,\star }^{f,V})\), where the two conserving schemes show excellent volume conservation properties in practice, with the observed relative volume changes always being less than \(0.1\%\). Finally, for general curvature flows of the form (1.6), we propose the linear scheme \(({\mathscr {A}}_m^F)^h\), which asymptotically exhibits equidistributed mesh points.
References
Altschuler, S.J., Wu, L.F.: Translating surfaces of the non-parametric mean curvature flow with prescribed contact angle. Calc. Var. Partial Differ. Equ. 2, 101–111 (1994)
Andrews, B.: Gauss curvature flow: the fate of the rolling stones. Invent. Math. 138, 151–161 (1999)
Andrews, B.: Motion of hypersurfaces by Gauss curvature. Pac. J. Math. 195, 1–34 (2000)
Athanassenas, M.: Volume-preserving mean curvature flow of rotationally symmetric surfaces. Comment. Math. Helv. 72, 52–66 (1997)
Barrett, J.W., Garcke, H., Nürnberg, R.: On the variational approximation of combined second and fourth order geometric evolution equations. SIAM J. Sci. Comput. 29, 1006–1041 (2007)
Barrett, J.W., Garcke, H., Nürnberg, R.: A parametric finite element method for fourth order geometric evolution equations. J. Comput. Phys. 222, 441–462 (2007)
Barrett, J.W., Garcke, H., Nürnberg, R.: On the parametric finite element approximation of evolving hypersurfaces in \({\mathbb{R}}^3\). J. Comput. Phys. 227, 4281–4307 (2008)
Barrett, J.W., Garcke, H., Nürnberg, R.: Parametric approximation of Willmore flow and related geometric evolution equations. SIAM J. Sci. Comput. 31, 225–253 (2008)
Barrett, J.W., Garcke, H., Nürnberg, R.: Finite element approximation of coupled surface and grain boundary motion with applications to thermal grooving and sintering. Eur. J. Appl. Math. 21, 519–556 (2010)
Barrett, J.W., Garcke, H., Nürnberg, R.: Variational Discretization of Axisymmetric Curvature Flows (2018). arXiv:1805.04322
Barrett, J.W., Garcke, H., Nürnberg, R.: Finite element methods for fourth order axisymmetric geometric evolution equations. J. Comput. Phys. 376, 733–766 (2019)
Basa, P., Schön, J.C., Salamon, P.: The use of Delaunay curves for the wetting of axisymmetric bodies. Q. Appl. Math. 52, 1–22 (1994)
Bernoff, A.J., Bertozzi, A.L., Witelski, T.P.: Axisymmetric surface diffusion: dynamics and stability of self-similar Pinchoff. J. Stat. Phys. 93, 725–776 (1998)
Cabezas-Rivas, E., Sinestrari, C.: Volume-preserving flow by powers of the \(m\)th mean curvature. Calc. Var. Partial Differ. Equ. 38, 441–469 (2010)
Coleman, B.D., Falk, R.S., Moakher, M.: Space-time finite element methods for surface diffusion with applications to the theory of the stability of cylinders. SIAM J. Sci. Comput. 17, 1434–1448 (1996)
Cox, G., Lowengrub, J.: The effect of spontaneous curvature on a two-phase vesicle. Nonlinearity 28, 773–793 (2015)
Deckelnick, K., Dziuk, G., Elliott, C.M.: Error analysis of a semidiscrete numerical scheme for diffusion in axially symmetric surfaces. SIAM J. Numer. Anal. 41, 2161–2179 (2003)
Deckelnick, K., Dziuk, G., Elliott, C.M.: Computation of geometric partial differential equations and mean curvature flow. Acta Numer. 14, 139–232 (2005)
Deckelnick, K., Schieweck, F.: Error analysis for the approximation of axisymmetric Willmore flow by \(C^1\)-finite elements. Interfaces Free Bound. 12, 551–574 (2010)
Dziuk, G.: An algorithm for evolutionary surfaces. Numer. Math. 58, 603–611 (1991)
Dziuk, G.: Convergence of a semi-discrete scheme for the curve shortening flow. Math. Models Methods Appl. Sci. 4, 589–606 (1994)
Dziuk, G., Kawohl, B.: On rotationally symmetric mean curvature flow. J. Differ. Equ. 93, 142–149 (1991)
Elliott, C.M., Fritz, H.: On approximations of the curve shortening flow and of the mean curvature flow based on the DeTurck trick. IMA J. Numer. Anal. 37, 543–603 (2017)
Finn, R.: Equilibrium Capillary Surfaces. Grundlehren der Mathematischen Wissenschaften, vol. 284. Springer, New York (1986)
Firey, W.J.: Shapes of worn stones. Mathematika 21, 1–11 (1974)
Ganesan, S., Tobiska, L.: An accurate finite element scheme with moving meshes for computing 3D-axisymmetric interface flows. Int. J. Numer. Methods Fluids 57, 119–138 (2008)
Gerhardt, C.: Flow of nonconvex hypersurfaces into spheres. J. Differ. Geom. 32, 299–314 (1990)
Grayson, M.A.: The heat equation shrinks embedded plane curves to round points. J. Differ. Geom. 26, 285–314 (1987)
Hartley, D.: Stability of near cylindrical stationary solutions to weighted-volume preserving curvature flows. J. Geom. Anal. 26, 2169–2203 (2016)
Hu, W.F., Kim, Y., Lai, M.C.: An immersed boundary method for simulating the dynamics of three-dimensional axisymmetric vesicles in Navier–Stokes flows. J. Comput. Phys. 257, 670–686 (2014)
Huisken, G.: The volume preserving mean curvature flow. J. Reine Angew. Math. 382, 35–48 (1987)
Huisken, G.: Asymptotic behavior for singularities of the mean curvature flow. J. Differ. Geom. 31, 285–299 (1990)
Huisken, G., Ilmanen, T.: The inverse mean curvature flow and the Riemannian Penrose inequality. J. Differ. Geom. 59, 353–437 (2001)
Jeffres, T.D.: Gauss curvature flow on surfaces of revolution. Adv. Geom. 9, 189–197 (2009)
Kovács, B., Li, B., Lubich, C.: A Convergent Evolving Finite Element Algorithm for Mean Curvature Flow of Closed Surfaces (2018). arXiv:1805.06667
LeCrone, J.: Stability and bifurcation of equilibria for the axisymmetric averaged mean curvature flow. Interfaces Free Bound. 16, 41–64 (2014)
Mantegazza, C.: Lecture notes on mean curvature flow. In: Progress in Mathematics, vol. 290. Birkhäuser/Springer Basel AG, Basel (2011)
Matioc, B.V.: Boundary value problems for rotationally symmetric mean curvature flows. Arch. Math. (Basel) 89, 365–372 (2007)
Mayer, U.F., Simonett, G.: A numerical scheme for axisymmetric solutions of curvature-driven free boundary problems, with applications to the Willmore flow. Interfaces Free Bound. 4, 89–109 (2002)
McCoy, J.A., Mofarreh, F.Y.Y., Wheeler, V.M.: Fully nonlinear curvature flow of axially symmetric hypersurfaces. NoDEA Nonlinear Differ. Equ. Appl. 22, 325–343 (2015)
McCoy, J.A., Mofarreh, F.Y.Y., Williams, G.H.: Fully nonlinear curvature flow of axially symmetric hypersurfaces with boundary conditions. Ann. Mat. Pura Appl. 4(193), 1443–1455 (2014)
Mikula, K., Ševčovič, D.: Evolution of plane curves driven by a nonlinear function of curvature and anisotropy. SIAM J. Appl. Math. 61, 1473–1501 (2001)
Urbas, J.I.E.: On the expansion of starshaped hypersurfaces by symmetric functions of their principal curvatures. Math. Z. 205, 355–372 (1990)
Veerapaneni, S.K., Gueyffier, D., Biros, G., Zorin, D.: A numerical method for simulating the dynamics of 3D axisymmetric vesicles suspended in viscous flows. J. Comput. Phys. 228, 7233–7249 (2009)
Zhao, Q.: A Sharp-Interface Model and Its Numerical Approximation for Solid-State Dewetting with Axisymmetric Geometry (2017). arXiv:1711.02402
Acknowledgements
The authors gratefully acknowledge the support of the Regensburger Universitätsstiftung Hans Vielberth.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
A Derivation of (2.17)
Here we demonstrate that (2.25b) and (2.26b) weakly impose (2.17). First we consider (2.25b) and the case \(\rho _0 = 0 \in \partial _0 I\).
We assume for almost all \(t\in (0,T)\) that \(\vec {x}(t) \in [C^1({\overline{I}})]^2\) and \(\varkappa _S(t) \in L^\infty (I)\). These assumptions and (2.3) imply that
for \({\overline{\rho }}\) sufficiently small, and for almost all \(t\in (0,T)\).
Let \(t \in (0,T)\). For a fixed \({\overline{\rho }} > 0\) and \(\varepsilon \in (0,{\overline{\rho }})\), we define
It follows from (A.1) that \((\vec {x}\cdot \,\vec {e}_1)\,\vec {\eta }_\varepsilon \) is integrable in the limit \(\varepsilon \rightarrow 0\). On choosing \(\vec {\eta } = \vec {\eta }_\varepsilon \in \underline{V}_\partial \) in (2.25b), we obtain in the limit \(\varepsilon \rightarrow 0\) that
Applying Fubini’s theorem and noting (A.1), as well as the boundedness of \(|\vec {x}_\rho |\) and \(\varkappa _{{\mathscr {S}}}\), yields the existence of a constant M such that
On the other hand, the right hand side in (A.2) converges to \((\vec {x}_\rho (0,t)\cdot \,\vec {e}_2)\,|\vec {x}_\rho (0,t)|^{-1}\) as \({\overline{\rho }}\rightarrow 0\), on recalling the smoothness assumptions on \(\vec {x}\). Combining this with (A.3) and (2.3) yields the boundary condition (2.17) for \(\rho = 0 \in \partial _0 I\). The proof for \(\rho = 1 \in \partial _0 I\) is analogous. Finally, the proof for (2.26b) is easily adapted from the above, on assuming that \(\vec {\varkappa }_{{\mathscr {S}}}(t) \in [L^\infty (I)]^2\) for almost all \(t\in (0,T)\).
B Existence proof for \(({\mathscr {A}}_m^f)^h\) and \(({\mathscr {A}}_m^{f,V})^h\)
We adapt [7, (2.12)–(2.14)] to \(({\mathscr {A}}_m^f)^h\) and \(({\mathscr {A}}_m^{f,V})^h\).
Theorem B.1
Let \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfy the assumptions \(({\mathfrak {A}})\) and \(({\mathfrak {B}})^h\), and assume that \(f{:}\,(a,b)\rightarrow {{\mathbb {R}}}\) with \(-\infty \le a<0<b\le \infty \) is strictly monotonically increasing, continuous and such that \(f((a,b))={\mathbb {R}}\). If \(b = -a = \infty \), then there exists a solution \((\delta \vec {X}^{m+1},\kappa ^{m+1}) \in \underline{V}^h_\partial \times V^h\) to \(({\mathscr {A}}_m^f)^h\) and \(({\mathscr {A}}_m^{f,V})^h\). Moreover, for general \(a<b\) there exists at most one solution.
Proof
Let \(f^{-1}{:}\,{{\mathbb {R}}}\rightarrow (a,b)\) denote the inverse of f. It follows from (4.13) and (4.15), on recalling (4.5), that
where \(g_1^m \in {{\mathbb {R}}}\) and \(g_0^m \in W^h_{\partial _0} \subset V^h\) are independent of \(\delta \vec {X}^{m+1}\) and \(\kappa ^{m+1}\). In particular, \(g_1^m = 0\) for \(({\mathscr {A}}_m^f)^h\). Substituting (B.1) into (4.3b) yields
where \(\ell ^m{:}\,\underline{V}^h\rightarrow {{\mathbb {R}}}\) is a linear functional defined by
It follows that (B.2) is the Euler–Lagrange variation of the minimization problem \(\min _{\vec {\eta } \in \underline{V}^h_\partial } {\mathscr {J}}^h(\vec {\eta })\), where
where \(\varPhi \in C^1({{\mathbb {R}}})\) denotes an antiderivative of \(f^{-1}\). We note that \(\varPhi {:}\,{{\mathbb {R}}}\rightarrow {{\mathbb {R}}}\) is strictly convex with \(\varPhi '(f(0))=f^{-1}(f(0))=0\) and hence we obtain that \(\varPhi \) is bounded from below and is coercive.
In the following we establish that the continuous functional \({\mathscr {J}}^h{:}\,\underline{V}^h_\partial \rightarrow {{\mathbb {R}}}\) is coercive, i.e. that \({\mathscr {J}}^h(\vec {\eta }) \rightarrow \infty \) as \(\Vert \vec {\eta }\Vert \rightarrow \infty \), where \(\Vert \cdot \Vert \) is a fixed norm on \(\underline{V}^h\). The main task is to bound the growth of the linear term \(\ell ^m(\vec {\eta })\) in terms of the first two terms in (B.3).
If \(b = -a = \infty \), it is possible to show that for all \(N \in {{\mathbb {N}}}\) there exists a positive constant \(C_0(N)\) such that
To see this, for \(N\in {{\mathbb {N}}}\) choose an \(R\in {{\mathbb {R}}}\) such that \(\min \{f(R),-f(-R)\} \ge N\), and define \(\Vert f\Vert _{\infty ,R} := \max _{q\in [-R,R]} |f(q)|\). On assuming without loss of generality that \(\varPhi (r) = \int _0^r f(q)\; {\mathrm{d}}q\), it holds for \(r > R\) that
In addition, for \(r \in [0,R]\) it trivially holds that
Combining (B.5) and (B.6) yields (B.4) for \(r \ge 0\). The case \(r\le 0\) can be treated analogously.
Given \(\vec {\eta } \in \underline{V}^h_\partial \), we define \(\vec {\zeta } = \vec {\eta } + \vec {{\mathfrak {f}}}^m \in \underline{V}^h\) with \(\vec {{\mathfrak {f}}}^m = \varDelta t_m\,g_1^m\,\vec {\pi }^h[ |\vec {\omega }^m|^{-2}\,\vec {\omega }^m] \in \underline{V}^h\). Then it holds for every \(N\in {{\mathbb {N}}}\) that
where, here and throughout, constants of the form \(C_i\) are independent of \(\vec {\zeta }\), but may depend on the data \(\vec {X}^m\), \(\vec {{\mathfrak {f}}}^m\) etc. Similarly, constants of the form \(C_i(N)\) may also depend on N, recall (B.4), but are independent of \(\vec {\zeta }\). On defining \(\int \!\!\!\!\!\mathrm{-}\eta = \frac{(\eta , |\vec {X}^m_\rho |)}{(1, |\vec {X}^m_\rho |)}\), and extending the definition to vector valued functions, it follows from (B.7) that
It remains to bound \(-{ III}-{ IV}-V\) from below. We have from the assumption \(({\mathfrak {B}})^h\) that \({ II} - V \ge N\,C_3\,|\textstyle \int \!\!\!\!\!\mathrm{-}\displaystyle \vec {\zeta }|\)\(- C_4\,|\textstyle \int \!\!\!\!\!\mathrm{-}\displaystyle \vec {\zeta }|\). Choosing \(N \ge 2\,C_4 / C_3\) implies that
In addition, it holds that
Choosing \((N+1)\,\delta \le \tfrac{1}{8}\) in (B.10) and combining with (B.8) and (B.9) implies that
which proves the coercivity of \({\mathscr {J}}^h(\vec {\eta })\).
We now consider the uniqueness of a solution to (B.2). Let \(\delta \vec {X}^{(i)} \in \underline{V}^h_\partial \), \(i=1,2\) be two solutions to (B.2). Then they satisfy
As \(f^{-1}\) is strictly monotonically increasing it immediately follows from (B.11) that \(\delta \vec {X}^{(1)} -\delta \vec {X}^{(2)} = \vec {X}^c \in {{\mathbb {R}}}^2\), and hence, on recalling (3.2), that
Now the strict monotonicity of \(f^{-1}\) implies that \(\vec {X}^c\cdot \,\vec {\omega }^m(q_j) = 0\) for all \(j=0,\ldots ,J\), and so the assumption \(({\mathfrak {B}})^h\) yields that \(\vec {X}^c=\vec {0}\). This shows the uniqueness of a solution to \(({\mathscr {A}}_m^f)^h\) and \(({\mathscr {A}}_m^{f,V})^h\).
Theorem B.1 yields existence of a unique solution for the schemes \(({\mathscr {A}}_m^f)^h\) and \(({\mathscr {A}}_m^{f,V})^h\) in the case (1.4a). For the case (1.4b) we only obtain uniqueness of a solution.
C Additional fully discrete schemes
Apart from the fully discrete schemes discussed in detail in Sect. 4, in this appendix we state some fully discrete variants of schemes introduced in Sect. 3 that are used for numerical computations in Sect. 5. We omit a detailed analysis of these schemes, but the interested reader can find some analysis of these schemes in the preprint version of this manuscript, see [10].
We begin with a fully discrete approximation of \(({\mathscr {B}}_h)^h\), (3.16a, b).
\(({\mathscr {B}}_m)^h\) Let \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \((\delta \vec {X}^{m+1}, \vec {\kappa }^{m+1}) \in \underline{V}^h_\partial \times \underline{V}^h\), where \(\vec {X}^{m+1} = \vec {X}^m + \delta \vec {X}^{m+1}\), such that
where \(\vec {{\mathfrak {K}}}^{m}(\vec {\kappa }^{m+1}) \in \underline{V}^h\) is such that
Lemma C.1
Let \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfy the assumption \((\mathfrak A)\). There exists a unique solution \((\delta \vec {X}^{m+1}, \vec {\kappa }^{m+1}) \in \underline{V}^h_\partial \times \underline{V}^h\) to \(({\mathscr {B}}_m)^h\).
Proof
See [10]. \(\square \)
Linear discrete analogues of the schemes \(({\mathscr {C}}_h)^{(h)}\), (3.18a, b), and \(({\mathscr {D}}_h)^{(h)}\), (3.19a, b), are given as follows.
\(({\mathscr {C}}_m)^{(h)}\) Let \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \((\delta \vec {X}^{m+1}, \kappa _{{\mathscr {S}}}^{m+1}) \in \underline{V}^h_\partial \times W^h_{(\partial _0)}\), where \(\vec {X}^{m+1} = \vec {X}^m + \delta \vec {X}^{m+1}\), such that
\(({\mathscr {D}}_m)^{(h)}\) Let \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \((\delta \vec {X}^{m+1},\vec {\kappa }_{{\mathscr {S}}}^{m+1})\in \underline{V}^h_\partial \times \underline{W}_{(\partial _0)}^h\), where \(\vec {X}^{m+1} = \vec {X}^m + \delta \vec {X}^{m+1}\), such that
We make the following assumption.
Lemma C.2
Let \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfy the assumptions \(({\mathfrak {A}})\) and \(({\mathfrak {C}}_{(\partial _0)})^{(h)}\). Then there exists a unique solution \((\delta \vec {X}^{m+1},\kappa _{{\mathscr {S}}}^{m+1}) \in \underline{V}^h_\partial \times W^h_{(\partial _0)}\) to \(({\mathscr {C}}_m)^{(h)}\).
Proof
See [10]. \(\square \)
Lemma C.3
Let \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfy the assumption \((\mathfrak A)\). Then there exists a unique solution \((\delta \vec {X}^{m+1},\vec {\kappa }_{{\mathscr {S}}}^{m+1}) \in \underline{V}^h_\partial \times \underline{W}_{(\partial _0)}^h\) to \(({\mathscr {D}}_m)^{(h)}\).
Proof
See [10]. \(\square \)
A nonlinear, and unconditionally stable, fully discrete variant of \(({\mathscr {D}}_h)^{(h)}\), (3.19a, b), is given as follows.
\(({\mathscr {D}}_{m,\star })^{(h)}\) Let \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \((\delta \vec {X}^{m+1},\vec {\kappa }_{{\mathscr {S}}}^{m+1})\in \underline{V}^h_\partial \times \underline{W}_{(\partial _0)}^h\), where \(\vec {X}^{m+1} = \vec {X}^m + \delta \vec {X}^{m+1}\), such that
Theorem C.1
Let \(\partial _0 I = \emptyset \) and let \(\vec {X}^m \in \underline{V}^h\) satisfy the assumption \(({\mathfrak {A}})\). Then there exists a solution \((\delta \vec {X}^{m+1},\vec {\kappa }_{{\mathscr {S}}}^{m+1}) \in \underline{V}^h_\partial \times \underline{W}_{(\partial _0)}^h\) to \(({\mathscr {D}}_{m,\star })^{(h)}\), (C.4a, b), if \(\varDelta t_m < 3\, \min _{{\overline{I}}} (\vec {X}^m\cdot \,\vec {e}_1)^2\).
Proof
See [10]. \(\square \)
Theorem C.2
Let \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfy the assumption \((\mathfrak A)\). Let \((\vec {X}^{m+1},\vec {\kappa }_{{\mathscr {S}}}^{m+1})\) be a solution to \(({\mathscr {D}}_{m,\star })^{(h)}\). Then it holds that
Proof
The proof is analogous to the proof of Theorem 4.1. \(\square \)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Barrett, J.W., Garcke, H. & Nürnberg, R. Variational discretization of axisymmetric curvature flows. Numer. Math. 141, 791–837 (2019). https://doi.org/10.1007/s00211-018-1013-z
Received:
Revised:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00211-018-1013-z