We consider numerical approximation to the evolution of a surface , determined by a flow map , driven by both mean curvature flow and a reaction–diffusion problem on the surface. That is, the surface velocity is given by the velocity law
where is the mean curvature and is the normal vector of the surface , and where is a given function of the solution of the nonlinear reaction–diffusion equation
Numerical approximation to the mean curvature flow — i.e. the case in (1.1) — was first addressed in 1990 by Dziuk , who proposed a finite element method based on an approximation of the mean curvature flow as a heat equation on a surface, in which the moving nodes of the finite element mesh determine the approximate evolving surface. However, proving convergence of Dziuk’s method or of other surface finite element based methods, such as the method proposed by Barrett, Garcke & Nürnberg  based on different variational formulations, or the method of Elliott & Fritz  based on DeTurck’s trick of re-parametrizing the surface, has remained an open problem for the mean curvature flow of closed two-dimensional surfaces.
In  we proved the first optimal-order convergence estimates for semi- and full discretisations of mean curvature flow of closed surfaces, by considering the coupled system for the velocity law together with evolution equations for the geometric quantities and , the normal vector and mean curvature.
Many practical applications concern mean curvature flow coupled with surface partial differential equations (PDEs), for example: tumour growth [7, 8, 6, 2]; surface dissolution [19, 17] (also see [15, Section 10.4]); diffusion induced grain boundary motion [20, 10, 29]; sharp interface limit of a perturbed Allen–Cahn equation ; etc. Numerical approximation to such forced mean curvature flow coupled with surface PDEs have been considered in some of these papers. For curves, convergence of numerical methods for such coupled problems was proved for forced curve shortening flow in [28, 3].
Convergence results, up to our knowledge, have not yet been proved for forced mean curvature flow of closed surfaces without regularization. For a regularized version of forced mean curvature flow of closed surfaces optimal-order convergence results for semi- and full discretisations were obtained in  and , respectively.
In this paper, we extend the approach and techniques of our previous paper  to the forced mean curvature flow problem (1.1)–(1.2) as a coupled problem, together with evolution equations for the normal vector and mean curvature. These evolution equations, as compared to mean curvature flow, cf. , contain additional forcing terms depending on . We present two fully discrete evolving finite element algorithms for the obtained coupled system. The first algorithm combines the two terms in spatial discretization by an idea of  treating conservation laws on an evolving surface. The second algorithm discretises the two terms separately in spatial discretization by using the velocity law for and the approach in . Both algorithms use evolving surface finite elements for spatial discretization and linearly implicit backward difference formulae for time integration, and for both algorithms the moving nodes of a finite element mesh determine the approximate evolving surface.
We show that the second algorithm can be written in the same matrix-vector form as the method proposed in  for the mean curvature flow. Thus the convergence analysis in  applies directly to the present algorithm for forced mean curvature flow as well. This yields uniform in time, optimal-order norm convergence results for the semi- and full discretisations of forced mean curvature flow, using at least quadratic evolving surface FEM and linearly implicit backward difference formulae of order two to five. For the first algorithm, we show a similar optimal-order convergence estimate of the evolving surface finite element semi-discretisation.
Finally, we present numerical experiments to support and complement the theoretical results. We present convergence tests for both algorithms, and also present an experiment with the numerical simulation for tumour growth model, using the parameters in  for the sake of easy comparison.
2 Evolution equations for mean curvature flow driven by diffusion on the surface
2.1 Basic notions and notation
We consider the evolving two-dimensional closed surface as the image
of a smooth mapping such that is an embedding for every . Here, is a smooth closed initial surface, and . In view of the subsequent numerical discretization, it is convenient to think of as the position at time of a moving particle with label , and of as a collection of such particles.
The velocity at a point equals
For a known velocity field , the position at time of the particle with label
is obtained by solving the ordinary differential equation (2.1) from to for a fixed .
For a function (, ) we denote the material derivative (with respect to the parametrization ) as
On any regular surface , we denote by the tangential gradient of a function , and in the case of a vector-valued function , we let . We thus use the convention that the gradient of has the gradient of the components as column vectors. We denote by the surface divergence of a vector field on , and by the Laplace–Beltrami operator applied to ; see the review  or [16, Appendix A] or any textbook on differential geometry for these notions.
We denote the unit outer normal vector field to by . Its surface gradient contains the (extrinsic) curvature data of the surface . At every , the matrix of the extended Weingarten map,
is a symmetric matrix (see, e.g., [30, Proposition 20]
). Apart from the eigenvalue
with eigenvector, its other two eigenvalues are the principal curvatures and . They determine the fundamental quantities
where denotes the Frobenius norm of the matrix . Here, is called the mean curvature (as in most of the literature, we do not put a factor 1/2).
2.2 Evolution equations for normal vector and mean curvature of a surface moving under mean curvature flow driven by diffusion on the surface
Mean curvature flow driven by diffusion on the surface, or for short forced mean curvature flow, sets the velocity (2.1) of the surface to
where denotes the normal velocity, and where the term involving is an external forcing, with a given function . Here, is the solution of the non-linear reaction–diffusion equation on the surface :
where is a given function.
The geometric quantities and on the right-hand side of (2.3) satisfy the following evolution equations, which are modifications of the classical evolution equations for pure mean curvature flow (i.e., ), cf. Huisken .
For a regular surface moving under forced mean curvature flow, the normal vector and the mean curvature satisfy
which, in combination with from (2.3), then gives the stated evolution equation for .
2.3 The system of equations used for discretization
Similarly as in , collecting the above equations, we have reformulated forced mean curvature flow as the system of semi-linear parabolic equations (2.5)–(2.6) on the surface coupled to the velocity law (2.3) and the surface PDE (2.4). The numerical discretization is based on a weak formulation of (2.3)–(2.6), together with the velocity equation (2.1):
for all test functions and , , and with . Here, we use the Sobolev space . This system is complemented with the initial data , , and .
for (not requiring ). Using the geometric quantities in the velocity law (2.3), using the following expression for the divergence
Throughout the paper both the usual Euclidean scalar product for vectors and the Frobenius inner product for matrices (which equals to the Euclidean product using an arbitrary vectorisation) are denoted by a dot.
3 Evolving finite element semi-discretization
3.1 Evolving surface finite elements
We formulate the evolving surface finite element (ESFEM) discretization for the velocity law coupled with evolution equations on the evolving surface, following the description in [24, 23], which is based on  and . We use simplicial finite elements and continuous piecewise polynomial basis functions of degree , as defined in [11, Section 2.5].
We triangulate the given smooth initial surface by an admissible family of triangulations of decreasing maximal element diameter ; see  for the notion of an admissible triangulation, which includes quasi-uniformity and shape regularity. For a momentarily fixed , we denote by the vector in that collects all nodes
of the initial triangulation. By piecewise polynomial interpolation of degree, the nodal vector defines an approximate surface that interpolates in the nodes . We will evolve the th node in time, denoted with , and collect the nodes at time in a column vector
We just write for when the dependence on is not important.
By piecewise polynomial interpolation on the plane reference triangle that corresponds to every curved triangle of the triangulation, the nodal vector defines a closed surface denoted by . We can then define globally continuous finite element basis functions
which have the property that on every triangle their pullback to the reference triangle is polynomial of degree , and which satisfy at the nodes for all These functions span the finite element space on ,
For a finite element function , the tangential gradient is defined piecewise on each element.
The discrete surface at time is parametrized by the initial discrete surface via the map defined by
which has the properties that for , that for all , and
The discrete velocity at a point is given by
In view of the transport property of the basis functions , the discrete velocity equals, for ,
where the dot denotes the time derivative . Hence, the discrete velocity is in the finite element space , with nodal vector .
The discrete material derivative of a finite element function with nodal values is
3.2 ESFEM spatial semi-discretizationss
Now we will describe the semi-discretization of the coupled system using both formulations of the surface PDE.
The finite element spatial semi-discretization of the weak coupled parabolic system (2.7) and (2.8) reads as follows: Find the unknown nodal vector and the unknown finite element functions and , , and such that, by denoting and ,
for all , , , and , with the surface given by the differential equation
The initial values for the nodal vector are taken as the positions of the nodes of the triangulation of the given initial surface . The initial data for , and are determined by Lagrange interpolation of , and , respectively.
Alternatively, the finite element spatial semi-discretization of the weak coupled parabolic system (2.7) and (2.9) determines the same unknown functions, but, instead of (3.2), the equations (3.1) and the ODE (3.3) are coupled to
for all .
In the above approaches, the discretization of the evolution equations for , and is done in the usual way of evolving surface finite elements. The velocity law (2.3) is enforced by a Ritz projection to the finite element space on . Taking instead just the finite element interpolation or projection does not appear to be sufficient for a convergence analysis (but in our experience both work well in practice). Note that the finite element functions and are not the normal vector and the mean curvature of the discrete surface .
3.3 Matrix–vector formulations
We collect the nodal values in column vectors , , and .
We define the surface-dependent mass matrix and stiffness matrix on the surface determined by the nodal vector :
with the finite element nodal basis functions . We define the nonlinear functions , and , where
with and given by, with the notations and ,
for and .
We further let, for (with the identity matrices )
When no confusion can arise, we will write for , for , and for .
can be written in the following matrix–vector form:
can be written in the matrix–vector formulation, where with same as before and with being the vector corresponding to the right-hand side of (3.4),
As in  and [23, Section 3.4], we compare functions on the exact surface with functions on the discrete surface , via functions on the interpolated surface , where denotes the nodal vector collecting the grid points on the exact surface, where are the nodes of the discrete initial triangulation .
Any finite element function on the discrete surface, with nodal values , is associated with a finite element function on the interpolated surface with the exact same nodal values. This can be further lifted to a function on the exact surface by using the lift operator , mapping a function on the interpolated surface to a function on the exact surface , provided that they are sufficiently close, see [12, 11].
Then the composed lift maps finite element functions on the discrete surface to functions on the exact surface via the interpolated surface is denoted by
4 Convergence of the semi-discretizations
We are now in the position to formulate the first main result of this paper, which yields optimal-order error bounds for the finite element semi-discretization (using finite elements of polynomial degree ) (3.1), and (3.2) or (3.4), with (3.3) of the system for forced mean curvature equations (2.7), and one of the weak formulations (2.8) or (2.9) for the surface PDE, with the ODE (2.1) for the positions. We introduce the notation
Consider the space discretizations (3.5) or (3.6) (in their matrix–vector form) of the coupled forced mean curvature flow problem (2.7), and (2.8) or (2.9), with (2.1), using evolving surface finite elements of polynomial degree . Suppose that the problem admits an exact solution that is sufficiently smooth on the time interval , and that the flow map is non-degenerate so that is a regular surface on the time interval .
Then for both formulations, there exists a constant such that for all mesh sizes the following error bounds for the lifts of the discrete position, velocity, normal vector and mean curvature hold over the exact surface for :
where the constant is independent of and , but depends on bounds of higher derivatives of the solution of the forced mean curvature flow and on the length of the time interval.
Sufficient regularity assumptions are the following: with bounds that are uniform in , we assume , , and for we have .
The proof of Theorem 4.1 for the semi-discretisation (3.6) follows from the stability bound [23, Proposition 7.1], combined with consistency of the numerical scheme, which are shown exactly as in [23, Lemma 8.1].
For the semi-discretisation (3.5) a stability proof can be obtained by combining the stability proofs of our previous works  and . The stability of the scheme is obtained by combining the results of [24, Proposition 6.1] (in particular part (A)) for the surface PDE, and of [23, Proposition 7.1] for the velocity law and for the geometric quantities. Numerical experiments presented in Section 7 demonstrate that semi-discrete error estimates, similar to those in Theorem LABEL:theorem:_coupled_error_estimate, also hold for the scheme (5.3).
5 Linearly implicit full discretization
For the time discretization of the system of ordinary differential equations of Section 3.3 we use a -step linearly implicit backward difference formula (BDF) with . For a step size , and with , let us introduce, for ,
|the discrete time derivative||(5.1)|
|the extrapolated value||(5.2)|
where the coefficients are given by and , respectively.
For (3.5) we obtain