# A convergent algorithm for mean curvature flow driven by diffusion on the surface

The evolution of a closed two-dimensional surface driven by both mean curvature flow and a reaction–diffusion process on the surface is formulated into a system, which couples the velocity law not only to the surface partial differential equation but also to the evolution equations for the geometric quantities, namely the normal vector and the mean curvature on the surface. Two algorithms are considered for the obtained system. Both methods combine surface finite elements as a space discretisation and linearly implicit backward difference formulae for time integration. Based on our recent results for mean curvature flow, one of the algorithms directly admits a convergence proof for its full discretisation in the case of finite elements of polynomial degree at least two and backward difference formulae of orders two to five. Numerical examples are provided to support and complement the theoretical convergence results (demonstrating the convergence properties of the method without error estimate), and demonstrate the effectiveness of the methods in simulating a three-dimensional tumour growth model.

## Authors

• 12 publications
• 15 publications
• 11 publications
07/22/2021

### A convergent finite element algorithm for mean curvature flow in higher codimension

Optimal-order uniform-in-time H^1-norm error estimates are given for sem...
10/21/2020

### A convergent finite element algorithm for generalized mean curvature flows of closed surfaces

An algorithm is proposed for generalized mean curvature flow of closed t...
04/10/2021

### Finite element error analysis for a system coupling surface evolution to diffusion on the surface

We consider a numerical scheme for the approximation of a system that co...
08/26/2019

### Some algorithms for the mean curvature flow under topological changes

This paper considers and proposes some algorithms to compute the mean cu...
01/19/2021

### Surfactant-dependent contact line dynamics and droplet adhesion on textured substrates: derivations and computations

We study the adhesion of a droplet with insoluble surfactant laid on its...
09/07/2019

10/19/2020

### Mean curvature motion of point cloud varifolds

This paper investigates a discretization scheme for mean curvature motio...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

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

 v=−Hν+g(w)ν, (1.1)

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

 ∂∙w+w∇Γ[X]⋅v−ΔΓ[X]w=F(w,∇Γ[X]w) on Γ[X], (1.2)

with a given nonlinear function . Problem (1.1)–(1.2) can be viewed as forced mean curvature flow.

Numerical approximation to the mean curvature flow — i.e. the case in (1.1) — was first addressed in 1990 by Dziuk [13], 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 [4] based on different variational formulations, or the method of Elliott & Fritz [18] 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 [23] 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 [27]; 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 [24] and [25], respectively.

In this paper, we extend the approach and techniques of our previous paper [23] 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. [22], 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 [14] 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 [23]. 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 [23] for the mean curvature flow. Thus the convergence analysis in [23] 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 [2] 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

 Γ[X]=Γ[X(⋅,t)]={X(p,t):p∈Γ0},

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

 ∂tX(p,t)=v(X(p,t),t). (2.1)

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

 ∂∙u(x,t)=ddtu(X(p,t),t) % for  x=X(p,t).

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 [9] 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,

 A(x)=∇Γν(x),

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

 H:=tr(A)=κ1+κ2,|A|2=κ21+κ22, (2.2)

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

 v=−Hν+g(w)ν=−Vν, (2.3)

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 :

 ∂∙w+w∇Γ[X]⋅v−ΔΓ[X]w=F(w,∇Γ[X]w), on Γ[X], (2.4)

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 [22].

###### Lemma 2.1.

For a regular surface moving under forced mean curvature flow, the normal vector and the mean curvature satisfy

 ∂∙ν =ΔΓ[X]ν+|A|2ν−∇Γ[X](g(w)), (2.5) ∂∙H =ΔΓ[X]H+|A|2H−ΔΓ[X](g(w))−|A|2g(w). (2.6)
###### Proof.

By using the normal velocity in the proof of [22, Lemma 3.3], or see also [5, Lemma 2.37], the following evolution equation for the normal vector holds:

 ∂∙ν=∇Γ[X]V.

On any surface , it holds true that (see [16, (A.9)] or [30, Proposition 24])

 ∇Γ[X]H=ΔΓ[X]ν+|A|2ν,

which, in combination with from (2.3), then gives the stated evolution equation for .

By revising the proof of [22, Theorem 3.4 and Corollary 3.5], or see [5, Lemma 2.39], with the normal velocity we obtain

 ∂∙H=ΔΓ[X]V+|A|2V,

which, again combined with from (2.3), yields the evolution equation for . ∎

### 2.3 The system of equations used for discretization

Similarly as in [23], 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):

 ∫Γ[X]∇Γ[X]v⋅∇Γ[X]φv+∫Γ[X]v⋅φv=−∫Γ[X]∇Γ[X]((H−g(w))ν)⋅∇Γ[X]φv ∫Γ[X]∇Γ[X]v⋅∇Γ[X]φv+∫Γ[X]v⋅φv=−∫Γ[X](H−g(w))ν⋅φv (2.7a) ∫Γ[X]∂∙ν⋅φν+∫Γ[X]∇Γ[X]ν⋅∇Γ[X]φν=∫Γ[X]|∇Γ[X]ν|2ν⋅φν−∫Γ[X]∇Γ[X](g(w))⋅φν (2.7b) ∫Γ[X]∂∙HφH+∫Γ[X]∇Γ[X]H⋅∇Γ[X]φH=∫Γ[X]|∇Γ[X]ν|2(H−g(w))φH ∫Γ[X]∂∙HφH+∫Γ[X]∇Γ[X]H⋅∇Γ[X]φH=+∫Γ[X]∇Γ[X](g(w))⋅∇Γ[X]φH, (2.7c)
 ddt∫Γ[X]wφw+∫Γ[X]∇Γ[X]w⋅∇Γ[X]φw=∫Γ[X]F(w,∇Γ[X]w)φw, (2.8)

for all test functions and , , and with . Here, we use the Sobolev space . This system is complemented with the initial data , , and .

An alternative weak formulation of (2.8), which is similar to (2.7b)–(2.7), but it involves the surface velocity, is based on

 ∫Γ[X]∂∙wφw+∫Γ[X]∇Γ[X]w⋅∇Γ[X]φw=∫Γ[X]F(w,∇Γ[X]w)φw−∫Γ[X](∇Γ[X]⋅v)wφw,

for (not requiring ). Using the geometric quantities in the velocity law (2.3), using the following expression for the divergence

 ∇Γ[X]⋅v= ∇Γ[X]⋅(−Vν)=−(∇Γ[X]V)⋅ν−V(∇Γ[X]⋅ν) = −(∇Γ[X](H−g(w))⋅ν−(H−g(w))(∇Γ[X]⋅ν),

is further rewritten (into a similar form as (2.7b) and (2.7)) as

 ∫Γ[X]∂∙wφw+∫Γ[X]∇Γ[X]w⋅∇Γ[X]φw=∫Γ[X]F2(H,ν,w,∇Γ[X]H,∇Γ[X]ν,∇Γ[X]w)φw. (2.9)

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.

###### Remark 2.2.

We remark here, however tempting to include, forcing terms involving the tangential gradient of cannot be covered by the analysis presented below. The cited stability arguments would not work for the Laplacian term for in (2.6), or for the velocity law (2.7).

## 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 [12] and [11]. 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 [14] 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

 x(t)∈R3N.

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

 ϕi[x]:Γh[x]→R,i=1,…,N,

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 ,

 Sh[x]=Sh(Γh[x])=span{ϕ1[x],ϕ2[x],…,ϕN[x]}.

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

 Xh(ph,t)=N∑j=1xj(t)ϕj[x(0)](ph),ph∈Γ0h,

which has the properties that for , that for all , and

 Γh[x(t)]=Γ[Xh(⋅,t)]={Xh(ph,t):ph∈Γ0h}.

The discrete velocity at a point is given by

 ∂tXh(ph,t)=vh(Xh(ph,t),t).

In view of the transport property of the basis functions [14], the discrete velocity equals, for ,

 vh(x,t)=N∑j=1vj(t)ϕj[x(t)](x)% with  vj(t)=˙xj(t),

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

 ∂∙huh(x,t)=ddtuh(Xh(ph,t))=N∑j=1˙uj(t)ϕj[x(t)](x)atx=Xh(ph,t).

### 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 ,

 ∫Γh[x]∇Γh[x]vh⋅∇Γh[x]φvh+∫Γh[x]vh⋅φvh=−∫Γh[x]∇Γh[x](Vhνh)⋅∇Γh[x]φvh−∫Γh[x]Vhνh⋅φvh (3.1a) ∫Γh[x]∂∙hνh⋅φνh+∫Γh[x]∇Γh[x]νh⋅∇Γh[x]φνh=∫Γh[x]α2hνh⋅φνh−∫Γh[x]∇Γh[x](g(wh))⋅φνh (3.1b) ∫Γh[x]∂∙hHhφHh+∫Γh[x]∇Γh[x]Hh⋅∇Γh[x]φHh=∫Γh[x]α2hVhφHh+∫Γh[x]∇Γh[x](g(wh))⋅∇Γh[x]φHh, (3.1c)
 ddt∫Γh[x]whφwh+∫Γh[x]∇Γh[x]wh⋅∇Γh[x]φwh=∫Γh[x]F(wh,∇Γh[x]wh)φwh (3.2)

for all , , , and , with the surface given by the differential equation

 ∂tXh(ph,t)=vh(Xh(ph,t),t),ph∈Γ0h. (3.3)

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

 ∫Γh[x]∂∙hwhφwh+∫Γh[x]∇Γ[X]wh⋅∇Γh[x]φwh=∫Γh[x]F2(Hh,νh,wh,∇Γh[x]Hh,∇Γh[x]νh,∇Γh[x]wh)φwh, (3.4)

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 :

 M(x)|ij= ∫Γh[x]ϕi[x]ϕj[x],A(x)|ij= ∫Γh[x]∇Γh[x]ϕi[x]⋅∇Γh[x]ϕj[x],i,j=1,…,N,

with the finite element nodal basis functions . We define the nonlinear functions , and , where

with and given by, with the notations and ,

 f1(x,n,H,w)|j+(ℓ−1)N= ∫Γh[x]α2h(νh)ℓϕj[x]−∫Γh[x](∇Γh[x](g(wh)))ℓ⋅ϕj[x], f2(x,n,H,w)|j= ∫Γh[x]α2hVhϕj[x]+∫Γh[x]∇Γh[x](g(wh))⋅∇Γh[x]ϕj[x], g(x,n,H,w)|j+(ℓ−1)N= −∫Γh[x]Vh(νh)ℓϕj[x]−∫Γh[x]∇Γh[x](Vh(νh)ℓ)⋅∇Γh[x]ϕj[x], F(x,w)|j= ∫Γh[x]F(wh,∇Γh[x]wh)ϕj[x]

for and .

We further let, for (with the identity matrices )

 M[d](x)=Id⊗M(x),A[d](x)=Id⊗A(x),K[d](x)=Id⊗(M(x)+A(x)).

When no confusion can arise, we will write for , for , and for .

Then, the equations (3.1) and (3.2) with (3.3), denoting

 u=(nH)∈R4N,

can be written in the following matrix–vector form:

 K[3](x)v= g(x,u,w), (3.5) M[4](x)˙u+A[4](x)u= f(x,u,w), ddt(M(x)w)+A(x)w= F(x,w) ˙x= v

Similarly, the equations (3.1) and (3.4) with (3.3), denoting

 u=⎛⎜⎝nHw⎞⎟⎠∈R5N,

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),

 K[3](x)v= g(x,u), (3.6) M[5](x)˙u+A[5](x)u= F(x,u), ˙x= v.

The system (3.6) for forced mean curvature flow is formally the same as the matrix–vector form of the coupled system for non-forced mean curvature flow, derived in [23], cf. (3.4)–(3.5) therein.

### 3.4 Lifts

As in [24] 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

 wLh=(ˆwh)l.

## 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

 xLh(x,t)=XLh(p,t)∈Γh[x(t)]forx=X(p,t)∈Γ[X(⋅,t)].
###### Theorem 4.1.

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 :

 ∥xLh(⋅,t)−idΓ(t)∥H1(Γ(t))3≤ Chk, ∥vLh(⋅,t)−v(⋅,t)∥H1(Γ(t))3≤ Chk, ∥νLh(⋅,t)−ν(⋅,t)∥H1(Γ(t))3≤ Chk, ∥HLh(⋅,t)−H(⋅,t)∥H1(Γ(t))≤ Chk, ∥wLh(⋅,t)−w(⋅,t)∥H1(Γ(t))≤ Chk, and also ∥Xlh(⋅,t)−X(⋅,t)∥H1(Γ0)3≤ Chk,

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].

###### Remark 4.2.

For the semi-discretisation (3.5) a stability proof can be obtained by combining the stability proofs of our previous works [24] and [23]. 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 ˙un= 1τq∑j=0δjun−j,and (5.1) the extrapolated value ˜un= q−1∑j=0γjun−1−j, (5.2)

where the coefficients are given by and , respectively.

We determine the approximations to , to , to , and to (only if not already collected into ) by the linearly implicit BDF discretisation of both systems (3.5) and (3.6).

For (3.5) we obtain

 K(˜xn