Thibaut Métivet

Arnaud Sengers

Mourad Ismail

Mourad Ismaïl b Emmanuel Maitre

Thibaut Metivet
email: thibaut.metivet@inria.fr

Diffusion-redistanciation schemes for 2D and 3D constrained Willmore flow: application to the equilibrium shapes of vesicles

Keywords: Diffusion generated motion, High-order geometrical flow, Willmore flow, Volume and area preserving, Level set

come   L'archive ouverte pluridisciplinaire

Introduction

Many different modelling situations can be cast as mathematical problems where an interface motion is driven by the minimisation of a geometric energy under geometric constraints. This is the case of multiphase flows, image segmentation, or elastic interface modeling, to cite a few. The geometric quantities involved in energy and constraints are for instance mean or Gaussian curvature of interface, surface area or enclosed volume. The algorithms and methods developed in this article give an efficient framework to address such situations.

Our work was motivated by one of the modeling situations that we now describe. Vesicles are systems of two fluids separated by a bi-layer membrane of phospholipid molecules. These objects can be considered as a simple model for Red Blood Cells (RBC). Since the number of such molecules is constant, this kind of interface has constant area. Therefore, its shape is determined by high order energy, i.e. the mean curvature is minimised. Moreover, there is no exchange across this interface, so that the enclosed volume is constant. Mathematically, the problem of finding a surface minimising its mean curvature is the well known problem of Willmore [START_REF] Simon | Existence of surfaces minimizing the willmore functional[END_REF]. In this work however, we are interested in this minimisation under the conservation of area and enclosed volume.

The numerical simulation of vesicles involve the resolution of two-fluid flows (for the inner and outer fluids) and fluid-structure interactions (for the membrane-fluids interaction), which is quite challenging: as a sharp object, a singularity occurs across the membrane bringing stress jump, which should either be dealt with explicit jump conditions, or appropriate numerical spreading. Since the membrane energy is of high geometrical order, its gradient involves high order derivatives of the unknowns. In addition, the inextensibility of the immersed interface is usually accounted for using elastic tension energies with high modulus which makes the resulting numerical problem very stiff.

A lot of such numerical simulations have been carried out by different teams using many numerical methods. To mention some representative works, we can cite the dynamic molecular method [38], the boundary integral method [3], the phase field method (see for example [4,59,[START_REF] Du | Adaptive finite element method for a phase field bending elasticity model of vesicle membrane deformations[END_REF][START_REF] John S Lowengrub | Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fission[END_REF]6]) or level set method [START_REF] Biben | Tumbling of vesicles under shear flow within an advected-field approach[END_REF]35,37,[START_REF] Salac | A level set projection model of lipid vesicles in general flows[END_REF][START_REF] Doyeux | Simulation of two-fluid flows using a finite element/level set method. application to bubbles and vesicle dynamics[END_REF]. [26] also proposes a model for vesicles implemented using a Lattice Boltzmann Method.

In the context of finite-element methods coupled with level set technique we can cite [START_REF] Laadhari | Computing the dynamics of biomembranes by combining conservative level set and adaptive finite element methods[END_REF][START_REF] Doyeux | Simulation of two-fluid flows using a finite element/level set method. application to bubbles and vesicle dynamics[END_REF]. We can also mention the work [24] based on a finite-element method where the membrane is modelled as a necklace of small rigid particles.

In the work above, the nonlinear coupling between the fluid flow and the geometric description of the interface is usually made explicit, which leads to severe restriction on the time step during simulation. Or it could be solved implicitly by a Newton-type method, which increases dramatically the cost per iteration. In [START_REF] Cottet | A semi-implicit level set method for multiphase flows and fluid-structure interaction problems[END_REF], a semi-implicit scheme was proposed where an ad hoc diffusion equation was used as a predictor step for the future position of a drop of liquid or a simple elastic interface.

Our approach in this work has as primary aim to propose a systematic way to built a predictor of the position of an interface with constant area and enclosed volume, when it moves to minimise high geometrical order energies, such as its mean curvature. In order to build such a predictor, we extended the diffusionthresholding/redistancing approach introduced in [41] and extended in [START_REF] Grzhibovskis | A convolution thresholding scheme for the willmore flow[END_REF][START_REF] Esedoglu | Threshold dynamics for high order geometric motions[END_REF][START_REF] Esedoglu | Diffusion generated motion using signed distance functions[END_REF][START_REF] Kublik | Algorithms for area preserving flows[END_REF] to the case of area and enclosed volume conservation.

Outline. The paper is organised as follows: in section 2 we review diffusion-thresholding and diffusionredistancing schemes principles, and introduce a new methodology to write higher order motion schemes, such as the Willmore flow, that is more systematic than in [START_REF] Esedoglu | Diffusion generated motion using signed distance functions[END_REF]. Then we present a new efficient method to take into account the area and enclosed volume constraints. In section 3 we address the problem of volume and area conservation, introducing an explicit analytic method to project on the constraint set. Then a section is devoted to numerical aspects of the implementation within finite-element approaches. We investigate in detail how to choose an optimal time step and plot numerical convergence curves for a basic scheme and an enhanced one. To evaluate the performance of the diffusion-redistanciation scheme without rescaling, we study the convergence of a torus under unconstrained Willmore flow to the optimal Clifford torus. To conclude, we present some numerical illustrations of the computation of 2D and 3D equilibrium shapes for vesicles that match well with those obtained with classical numerical schemes.

Numerical schemes for higher geometrical motion of interfaces

The class of diffusion-generated motion was introduced by the work [41] of Merriman, Bence and Osher. They proposed an efficient algorithm for computing the mean curvature flow of a surface (a curve in two dimensions) without any direct computation of the mean curvature. The algorithm consists in repeating two steps, namely a diffusion step (also named convolution step) and a thresholding step. During the diffusion step, the characteristic function representing the surface is diffused for a certain time step, and the 1 2 iso-level moves proportionally to the local curvature. A characteristic function is then recovered by thresholding the resulting function at 1 2 , allowing to iterate the process.

Algorithm 1 Original Convolution-Thresholding scheme while t < t f do Solve ∂ t φ -∆φ = 0 with initial condition φ i = 1 (n) for time δt Construct the new characteristic function 1 (n+1) = 1 {φ≥ 1

The main advantage of this method is its simplicity and unconditional stability: the convolution step consists in the numerical resolution of a heat equation which can be achieved efficiently, while the standard phase field approach would involve solving a more complicated non-linear equation or computing the curvature and solving a transport equation.

Due to its conceptual simplicity and numerical efficiency, various extensions of the method have been proposed since its introduction. Extension to the multiphase case is trivial for symmetric junctions as one only has to diffuse each interface and threshold to the biggest one in each region. The case of nonsymmetric junctions and different surface tensions is treated in [START_REF] Steven | A diffusion-generated approach to multiphase motion[END_REF]. However the resulting algorithm relies on a spatially dependent thresholding that require essentially to locate triple junctions. Also the case where four or more phases meet is only treated empirically. In [START_REF] Esedoglu | Threshold dynamics for networks with arbitrary surface tensions[END_REF] a reformulation of the problem is done and the associated algorithm changes the diffusion step to take into account the difference between the surface tensions. An area preserving motion by mean curvature can be found in [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF] by changing the threshold value from 1 2 to a real number λ found using a Newton method. We refer to the next section for more details about area conservation.

Extensions of the convolution-thresholding method to the Willmore flow and other higher order geometric motion, such as surface diffusion, have been proposed independently in [START_REF] Esedoglu | Threshold dynamics for high order geometric motions[END_REF][START_REF] Grzhibovskis | A convolution thresholding scheme for the willmore flow[END_REF] in both two and three dimensions. We recall that the Willmore flow is given by the following normal velocity :

W =    ∆ Γ H + H 3 2 in 2D ∆ Γ H + 2H H 2 -K in 3D, (1)
where

H =    κ 1 in 2D κ 1 + κ 2 2 in 3D, and K = κ 1 in 2D κ 1 κ 2 in 3D,
are respectively the mean and Gaussian curvatures of the surface, κ i denoting the usual principal curvatures.

In 2D [START_REF] Esedoglu | Threshold dynamics for high order geometric motions[END_REF] and in 3D [START_REF] Grzhibovskis | A convolution thresholding scheme for the willmore flow[END_REF], a local expansion of the convoluted characteristic function is performed and shows that the velocity W appears at the second order of the diffused interface. More precisely, in the 2D case, any point of the interface can be relocated at the origin in a way that the normal to the curve is aligned with the y-axis. The behavior of the point of interest during the convolution step is then described by the following expansion:

F (0, y, δt) = 1 2 - 1 2 √ π yδt -1 2 + 1 2 √ π κδt 1 2 - 1 4 √ π W δt 3 2 + O δt 5 2 (2)
where F (0, y, δt) denotes the function obtained by diffusing the characteristic function for a time δt. Thresholding the result at F (0, y, δt) = 1 2 then results in a normal displacement of the point by κδt, which gives the classical mean curvature scheme. As the velocity term of interest for the Willmore flow can be found at the second order in expansion eq. (2), one can then compute this flow by extracting it with a linear combination of two different solutions of the convolution step taken at different time steps √ 2∆t/θ and θ √ 2∆t for θ = 1 that eliminates the curvature term. Upon thresholding at an appropriate level [START_REF] Esedoglu | Threshold dynamics for high order geometric motions[END_REF] found that the curve moved following the Willmore flow since the solution is

y = δtW + O δt 3 2 .
At the expense of an additional diffusion at each step and an error of order δt 3 2 , we can write an algorithm for the Willmore flow similar to the mean curvature one. A similar expansion in the 3D case can be found in [START_REF] Grzhibovskis | A convolution thresholding scheme for the willmore flow[END_REF] and yields the same expression to compute the Willmore flow.

This method seems highly valuable because it allows to compute a fourth-order flow without any differentiation, especially in the context of finite-element methods with low-order polynomial discretisations of the phase-field function φ.

Perform one step of ∂ t ψ -∆ψ = 0 with initial condition ψ 0 = 1 (n) for time steps

δt 1 = √ 2δt/θ and δt 2 = θ √ 2δt Compute D = (2δt) 1 4 [ √ θψ (δt 1 ) -1 √ θ ψ (δt 2 )]
Construct the new characteristic function

1 (n+1) = 1 {D≥ 1 2 (2δt) 1 4 √ θ-1 √ θ } t ← t + δt
However the convolution-thresholding method has some identified drawbacks. Its major issue is its inaccuracy and its propensity to get "stuck" when the mesh is not refined enough at the interface. More precisely, as pointed out in [START_REF] Merriman | Motion of multiple junctions: A level set approach[END_REF], given a fixed uniform grid, if the motion during an iteration is smaller than the grid resolution, the thresholding step can reset the interface to its initial configuration and loop indefinitely. A solution to this problem can be found in [48], where the authors use adaptive grids to refine the mesh near the interface where the resolution is more important. Such adaptive strategy is also essential for accuracy considerations, especially in the case of higher order flows, but can become algorithmically complicated due to its use of unequally spaced FFTs in the dynamic case or in dimension 3.

An alternative approach to the problem is to allow the function to contain subgrid informations. In [START_REF] Esedoglu | Diffusion generated motion using signed distance functions[END_REF], the authors replace to that end the characteristic function by a signed distance function and explore the motions that can be obtained through the diffusion and the redistanciation of a signed distance function. A similar expansion as eq. ( 2) of the solution of the heat equation starting from a distance function shows that the first order term is still a mean curvature term. At the expense of a redistanciation step, though more costly than the thresholding step, one can obtain significantly more accurate computations than the classical convolution thresholding algorithm. However, second order term in this case is not the desired Willmore term any more and an additional correction is therefore necessary to compute the Willmore flow.

From a theoretical point of view, convergence of the classical algorithm of motion by mean curvature is well understood and various proofs of the consistency of the scheme exist. Some are based on a comparison principle [START_REF] Lawrence | Convergence of an algorithm for mean curvature motion[END_REF]2,[START_REF] Chambolle | Convergence of an algorithm for the anisotropic and crystalline mean curvature flow[END_REF], other are based on a gradient flow interpretation [START_REF] Almgren | Curvature-driven flows: a variational approach[END_REF]36]. Introduced in the two papers [START_REF] Esedoglu | Threshold dynamics for networks with arbitrary surface tensions[END_REF][START_REF] Laux | Convergence of the thresholding scheme for multi-phase mean-curvature flow[END_REF], a new interpretation of the convolution-thresholding algorithm as a gradient minimising flow has initiated more recent works on the subject [31,[START_REF] Laux | Brakke's inequality for the thresholding scheme[END_REF], while [START_REF] Swartz | Convergence of diffusion generated motion to motion by mean curvature[END_REF] relies on a formal matched asymptotic method. The extension to the multiphase flow can be found in [START_REF] Laux | Convergence of the thresholding scheme for multi-phase mean-curvature flow[END_REF] and the case of volume-preserving motion is studied in [33,43]. The recent results on the gradient flow approach are reviewed in [START_REF] Laux | Gradient-flow techniques for the analysis of numerical schemes for multi-phase mean-curvature flow[END_REF]. However, the convergence of higher order motions, especially the Willmore flow, remains an open question.

In this section, we propose a new local expansion of the signed distance function near the interface. It yields similar results to the one in [START_REF] Esedoglu | Diffusion generated motion using signed distance functions[END_REF] for the 2D case but our approach seems more intrinsic and its generalisation to higher dimensions and orders is straightforward. We propose diffusion-redistanciation algorithms to compute the Willmore flow in 2D and 3D. As the term in the second order of the expansion is not quite the Willmore term, one has to add a correction term. The 3D version in particular shows a completely new approach and features a correction term using a diffusion of d 2 .

Expansion of a distance function convoluted by the heat kernel

As suggested by the previous work of Esedoglu, Ruuth and Tsai [START_REF] Esedoglu | Diffusion generated motion using signed distance functions[END_REF], the high-order geometrical motions of a surface can be computed using the convolution of a distance function to this interface by the heat equation kernel. In this section, we thus derive an expansion of the solution of the heat equation with a distance function as its initial condition. In contrast with [START_REF] Esedoglu | Diffusion generated motion using signed distance functions[END_REF], we compute this expansion for both the twoand three-dimensional cases using a fully implicit approach.

We consider d a signed distance function to some interface Γ, and introduce the heat equation:

∂ t φ -∆φ = 0. (3)
supplemented with the initial condition:

φ(t = 0) = d. (4)
The Taylor expansion of the heat equation solution with respect to time then reads:

φ(δt) = φ(0) + ∂ t φ(0) δt + ∂ tt φ(0) δt 2 2 + O δt 3 (5)
Recalling eq. (3) and the initial condition eq. (4), we then have ∂ t φ(0) = ∆d and ∂ tt φ(0) = ∆ 2 d. We denote w = ∆ 2 d 2 and can then rewrite eq. ( 5) as:

G δt ≡ φ(δt) = d + ∆d δt + ∆ 2 d 2 δt 2 + O δt 3 = d + ∆d δt + w δt 2 + O δt 3 . ( 6
)
Remark 1. When time discretisation of (3) will come into play, we will have to consider a second order scheme to get the right expansion. Using a first order Euler scheme, which could be desirable for its stability properties, leads to double the second order term magnitude. We refer to section 4.2 for more details.

By definition, ∆d is the total curvature of the local d iso-surface:

∆d = κ = (n -1)H
with n the space dimension and H the mean curvature. We can also express ∆ 2 d close to the interface Γ as a function of the interface geometrical quantities, and in particular the Laplace-Beltrami operator of the curvature ∆ Γ κ = (n -1)∆ Γ H. More specifically, if p denotes the projection of x on Γ, we have

∆ Γ f = ∆ (f (x -d∇d)) | Γ ≡ ∆ (f • p) | Γ for any f ∈ C 2 (Γ) [12]
, so that:

∆ Γ (∆d) = ∆(∆d • p)| Γ = ∇ i (∇ j p i ∇ j ∆d) | Γ = ∇ i ((δ ij -∇ i d ∇ j d -d∇ i ∇ j d) ∇ j ∆d) | Γ = -∆d ∇ N ∆d -∇ 2 N ∆d + ∆ 2 d -d∇ i (∇ i ∇ j d ∇ j ∆d) | Γ
where N ≡ ∇d is the interface normal and we have introduced the normal gradient ∇ N ≡ N • ∇. Applying the restriction to Γ ≡ {d = 0}, we finally obtain:

∆ 2 d = ∆ Γ κ + ∇ 2 N κ + κ ∇ N κ (7)
We can also compute:

∇ N κ = ∇ i d ∇ i ∇ k ∇ k d = ∇ i d ∇ k ∇ k ∇ i d = ∇ k ( ( ( ( ( ( ( ∇ i d ∇ k ∇ i d) -∇ k ∇ i d ∇ i ∇ k d = -Tr (∇ 2 d) 2 since the distance property |∇d| = 1 gives ∇ i d ∇ k ∇ i d = ∇ k |∇d| 2 /2 = 0. Recalling that ∇ 2 d is the
Weingarten map of Γ, we can then write:

∇ N κ = - i κ 2 i = -κ 2 2K -κ 2 = -H 2 in 2D -2 2H 2 -K in 3D (8)
where we have introduced the principal curvatures κ i , i = 1, . . . , n -1, and the Gauss curvature K = i κ i .

The second normal derivative can then be written as:

∇ 2 N κ = ∇ i d ∇ i ∇ N κ = -2 ∇ i d ∇ i ∇ j ∇ k d ∇ k ∇ j d = -2 [∇ j ( ( ( ( ( ( ( ∇ i d ∇ k ∇ i d ∇ j ∇ k d) -∇ j ∇ i d ∇ k ∇ i d ∇ j ∇ k d -( ( ( ( ( ( ∇ i d ∇ k ∇ i d ∇ j ∇ j ∇ k d] = 2 Tr (∇ 2 d) 3
that is:

∇ 2 N κ = 2 i κ 3 i = 2κ 3 2κ κ 2 -3K = 2H 3 in 2D 4H 4H 2 -3K in 3D. (9)
Inserting eq. ( 8) and eq. (9) into eq. ( 7), we finally obtain:

∆ 2 d ≡ 2w = ∆ Γ κ + κ 3 ∆ Γ κ + κ(κ 2 -4K) = ∆ Γ H + H 3 in 2D 2 ∆ Γ H + 4H(H 2 -K) in 3D. ( 10
)
As expected, we observe that the second-order term in the expansion of the heat equation solution (6) features the required high-order geometrical quantities for the geometrical flows we are interested in. In the following, we shall in particular devise specific algorithms for the simulation of Willmore flows, but the method presented here is generic, and can be carried on at higher orders similarly.

Convolution of a generic function of the distance

For practical reasons that will prove useful thereafter, we consider in this section the expansion of the heat equation solution eq. ( 3), but with the initial condition:

φ(t = 0) = f (d)
where f is some generic C ∞ (R) function. Similarly to eq. ( 6), we denote G δt the corresponding solution, and write as before:

G (f ) δt = f (d) + ∆f (d) δt + ∆ 2 f (d) δt 2 2 + O δt 3 . (11)
Recalling that |∇d| = 1, and using the usual derivation rules, we can then compute: 1) ∆d and:

∆f (d) = ∇ • (f (1) ∇d) = f (2) + f (
∆ 2 f (d) = f (4) (d) + 2f (3) (d)∆d + f (2) (d) (∆d) 2 + 2∇ N ∆d + f (1) (d)∆ 2 d = f (4) (d) + 2f (3) (d)κ + f (2) (d) κ 2 + 2∇ N κ + 2f (1) (d)w ≡ f (4) (d) + 2f (3) (d)κ + 2f (2) (d)c + 2f (1) (d)w
where we have introduced:

c ≡ κ 2 2 + ∇ N κ =        - κ 2 2 2K - κ 2 2 =    - H 2 2 in 2D -2(H 2 -K) in 3D. (12)
We eventually get:

G (f ) δt = f (d) + δt f (2) (d) + f (1) (d)κ + δt 2 f (4) (d) 2 + f (3) (d)κ + f (2) (d)c + f (1) (d)w + O δt 3 . ( 13
)

Expansion order selection

We now turn to the use of combinations of G for different time-steps in order to select the appropriate order in the expansion eq. (6). To minimise assembly and preconditioner computation costs, we restrict our analysis to the case of two different time-steps denoted αδt and βδt respectively. To simplify the notations, we denote in the following G α ≡ G αδt and G β ≡ G βδt . We then consider the general linear combination a G α + b G β which upon expansion reads:

a G α + b G β = (a + b) d + (aα + bβ) κδt + (aα 2 + bβ 2 ) wδt 2
The coefficients can then be adjusted to cancel the term of order δt 2 (the "Order 1" combination) or the term of order δt (the "Order 2" combination). As mentioned above, we want to use only α and β for both combinations; but the multiplicative coefficients a and b can be chosen arbitrarily. We thus introduce the "Order 1" (resp. "Order 2") coefficients a 1 and b 1 (resp. a 2 and b 2 ) and write the corresponding order selections as:

Order 1 : a 1 α + b 1 β = 1, a 1 α 2 + b 1 β 2 = 0 (14)
Order 2 :

a 2 α + b 2 β = 0, a 2 α 2 + b 2 β 2 = -1 (15)
Note that we have arbitrarily imposed a 1 α + b 1 β = 1 and a 2 α 2 + b 2 β 2 = -1 (these quantities actually only need to be non-zero) to directly recover the combinations used for our Willmore flow algorithms below.

With the additional constraint that α > 0, β > 0, the general solution is:

β = α √ -a 1 √ b 1 , b 2 = a 2 √ b 1 √ -a 1 .
We thus choose the following parameters:

α = √ 2, β = √ 2 2 a 1 = - √ 2 2 , b 1 = 2 √ 2 a 2 = -1, b 2 = 2
so that:

a 1 G α + b 1 G β = 3 √ 2 2 d + κδt + O δt 3 (16)
and:

a 2 G α + b 2 G β = d -wδt 2 + O δt 3 (17)

Willmore flow

In this section, we derive an algorithm which relies on solving heat equation problems starting from a distance function to compute high-order motions of an interface in a two-or three-dimensional domain. We focus our analysis more specifically to the Willmore flow, but the approach can be straightforwardly extended to a large class of surface diffusion motions.

We consider the Willmore energy functional [START_REF] Thomas | Riemannian geometry[END_REF]:

E W (Γ) ≡ Γ H 2 . ( 18
)
The surfaces which minimise eq. ( 18) over the set of one-or two-dimensional surfaces can be obtained from the Willmore flow gradient descent, which moves iteratively some initial surface with the normal velocity W N where (c.f. [START_REF] Thomas | Riemannian geometry[END_REF]):

W =    ∆ Γ H + H 3 2 in 2D ∆ Γ H + 2H(H 2 -K) in 3D.
For a surface Γ implicitly defined by the signed distance function as Γ = x ∈ R 3 , d( x) = 0 , the Willmore flow corresponds to a usual transport equation:

∂ t d + W N • ∇d = 0. ( 19
)
Recalling that ∇d = N , this amounts to solving ∂ t d+W = 0, or if we consider a first-order time-discretisation with time-step δt,

d (p+1) = d (p) -W δt. (20)
We can then observe that the right-hand side of eq. ( 20) can be obtained from the expansions of solutions to the heat equation derived above. More precisely, the second order term wδt 2 (eq. ( 10)) provides the required high-order surface term ∆ Γ H. This is reminiscent of the approach of [START_REF] Esedoglu | Threshold dynamics for high order geometric motions[END_REF][START_REF] Grzhibovskis | A convolution thresholding scheme for the willmore flow[END_REF], where convolutions of the characteristic function defining the interface are combined to evolve the interface according to the Willmore flow. However, while the second-order terms appearing in the expansion of the convolution of the characteristic function by the heat kernel provide the exact Willmore velocity, this is no more the case for the distance function, and we need to add some correction terms in order to retrieve the correct Willmore flow.

We also stress that d (n) -W δt is not in general a distance function, so that we need to add a redistanciation step to iterate the process. This redistanciation is the counterpart to the thresholding step in the usual convolution-thresholding algorithms for mean curvature [41] or higher-order [START_REF] Esedoglu | Threshold dynamics for high order geometric motions[END_REF][START_REF] Grzhibovskis | A convolution thresholding scheme for the willmore flow[END_REF] motions using a characteristic function. In the following, we denote Redist this redistanciation step, and will use a fast-marching algorithm to actually compute it. We provide more details about our implementation in section 4.

2D case. In 2D, the correction is rather simple, and only involves the curvature H 3 , namely:

W 2D = 2w - H 3 2 = 2 δt 2 (d -a 2 G α -b 2 G β ) - H 3 2 . ( 21
)
The curvature H can be computed with different methods, including using the order-1 selection eq. ( 16), the c correction term introduced above in eq. ( 12) or direct methods using ∆d. In our computations we will use either the order-1 diffusion combination method or the c correction term computed with diffusions of d 2 (c.f. eq. ( 23)) as for the three-dimensional case. We refer to section 4.3 for a discussion on the efficiency of each method.

As pointed out in remark 1 and discussed in more detail in section 4.2, the choice of the numerical scheme used to solve the heat equation has an impact on the development eq. ( 5). We present here algorithms using either an implicit Euler or a Crank-Nicolson scheme and choose our combination accordingly.

Our diffusion-redistanciation algorithm for the Willmore flow of a curve in two dimensions is finally given in alg. 3.

Algorithm 3 2D Willmore flow 1 while t < t f do Solve ∂ t φ -∆φ = 0 with initial condition φ i = d (n) for times t = √ 2δt and t = δt/ √ 2 → G √ 2δt , G δt/ √ 2 .
Compute the curvature H.

Compute D =        D E = -G √ 2δt + 2G δt/ √ 2 + H 3 2 δt 2 D CN = -G √ 2δt + 2G δt/ √ 2 + H 3 4 δt 2 . Compute the new signed distance function d (n+1) = Redist(D). t ←    t + δt 2 (Euler) t + δt 2 2 (Crank-Nicolson)
During our simulations, we however preferred to use an algorithm that computes the correction the same as in the 3D case method using the convolution of d 2 , presented in the next paragraph, as it proved more robust in some cases. The corresponding algorithm is described in alg. 4. Despite two additional convolutions, this algorithm seems to produce better results for our numerical experiments. It has the advantage of using a correction term which is already of order two. We suspect that multiplying the correction by δt 2 in the first version can have some detrimental numerical effects due to the space discretisation errors.

Algorithm 4 2D Willmore flow with alternative correction

1 while t < t f do Solve ∂ t φ -∆φ = 0 with initial condition φ i = d (n) for times t = √ 2δt and t = δt/ √ 2 → G √ 2δt , G δt/ √ 2 . Solve ∂ t φ -∆φ = 0 with initial condition φ i = d (n)2 for times t = √ 2δt and t = δt/ √ 2 → G d 2 √ 2δt , G d 2 δt/ √ 2 .
Compute the curvature H.

Compute D = Hd (n)2 4 + -G √ 2δt + 2G δt/ √ 2 • 1 -Hd (n) 2 + -G d 2 √ 2δt + 2G d 2 δt/ √ 2 • H 4 Compute the new signed distance function d (n+1) = Redist(D). t ←    t + δt 2 (Euler) t + δt 2 2 (Crank-Nicolson)
3D case. In 3D, the correction term is more complicated, as it involves both the mean and the Gaussian curvatures. More precisely, we have:

W 3D = w -2H(H 2 -K) = w + Hc
where c is the term introduced in eq. ( 12). To compute this correction term c, we note from eq. ( 13) that we need to choose f such that its second derivative does not vanish. To this end, we choose arbitrarily

f (d) = d 2 and introduce G d 2 α = d 2 + 2αδt (1 + dκ) + 2(αδt) 2 (c + dw) + O δt 3 . (22)
so that we can extract the correction term from HG d 2 α . However, using only G d 2 α to compute c requires to compensate the zero-and first-order terms in eq. ( 22), which can lead to large errors in particular regarding the 2δt d κ term. Note that while this latter term disappears as h → 0 on the interface (as d| Γ = 0 by definition), the discrete term brings up O (h) errors which can be large when in the δt order term. To prevent this, we choose to also select directly the secondorder term in eq. ( 22) using a linear combination as introduced in eq. [START_REF] Esedoglu | Diffusion generated motion using signed distance functions[END_REF]. We thus use:

a 2 G d 2 α + b 2 G d 2 β = d 2 -2δt 2 (c + dw) + O δt 3 . ( 23
)
Recalling W 3D = w + Hc, and combining eq. ( 17) and eq. ( 23), we finally obtain:

W 3D = 1 δt 2 d - Hd 2 2 -(a 2 G α + b 2 G β ) (1 -H d) - H 2 a 2 G d 2 α + b 2 G d 2 β ( 24
)
and our algorithm for the three-dimensional Willmore flow is given in alg. 5.

Remark 2. Note that the flow algorithms actually compute a gradient descent with step δt 2 . However, the evolution and diffusion time-steps can easily be decoupled involving the ratio

δt f low δt 2 dif f
, but our numerical tests suggest that keeping this ratio close to 1 leads to the more stable simulations. This decoupling can nevertheless prove useful to implement efficient backtracking methods for the choice of the flow time-step while avoiding unnecessary resolutions of the diffusion equation.

Diffusion-redistancing schemes with volume and area conservation

Since we are interested in computing shapes in R n , n = 2 or n = 3, we will in this chapter call volume the measure with respect to the Lebesgue measure of R n and area the surface measure of hypersurfaces of R n .

Algorithm 5 3D Willmore flow 1 while t < t f do Solve ∂ t φ -∆φ = 0 with initial condition φ i = d (n) for times t = √ 2δt and t = δt/ √ 2 → G √ 2δt , G δt/ √ 2 . Solve ∂ t φ -∆φ = 0 with initial condition φ i = d (n)2 for times t = √ 2δt and t = δt/ √ 2 → G d 2 √ 2δt , G d 2 δt/ √ 2 .
Compute the curvature H.

Compute D = Hd (n)2 2 + -G √ 2δt + 2G δt/ √ 2 • 1 -Hd (n) + -G d 2 √ 2δt + 2G d 2 δt/ √ 2 • H 2 Compute the new signed distance function d (n+1) = Redist(D). t ← t + 2δt 2 (Euler) t + δt 2 (Crank-Nicolson)

Review of existing methods

Problems involving area or volume conservation were addressed in the framework of diffusion-redistanciation schemes by [START_REF] Kublik | Algorithms for area preserving flows[END_REF], and their convergence properties in [33]. In [START_REF] Kublik | Algorithms for area preserving flows[END_REF], a new algorithm for area preserving flows in two dimensions is introduced by considering normal velocities

v N = κ -κ + S
where κ is the curvature, κ its average on the interface, and S is an additional term which depends on the application. The algorithm relies in 2D on a property linking the mean curvature with the surface area and the genus number, which in case of a connected hypersurface S of genus 0 enclosing an open bounded set Ω, boils down to

κ = 2π |S| .
Using the divergence theorem, one can then compute |S| with a volume integral as

|S| = Ω ∆d dx.
In the case where S = 0, to ensure an accurate volume conservation, a Newton method is applied to find a real number λ * close to κ computed above (which is used as initialisation) to correct the diffusion generated motion by raising or lowering the convolution of the signed distance function (as was introduced the framework of convolution-thresholding schemes by [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF]). In the case n = 3, while the trick to compute efficiently the mean curvature is not anymore valid, we could still apply the method by computing numerically the mean of the mean curvature.

The level-set community has also developed several techniques to address volume conservation [START_REF] Sussman | An improved level set method for incompressible two-phase flows[END_REF][START_REF] Sussman | An efficient, interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow[END_REF][START_REF] Russo | A remark on computing distance functions[END_REF] in the modelling of multiphase flows. This corresponds to correct the Hamilton-Jacobi equation of the redistanciation step by a term which ensures that conservation. Namely, the following equation:

∂ τ φ = sgn(φ 0 )(1 -∇φ ) + λf (φ)
is solved with λ computed so that Ωij H(φ)dx is conserved on each cell Ω ij of the grid, and f chosen to localise around the interface.

An even more straightforward approach is a post-or preprocessing trick due to Smolianski [55] in the level-set framework, where a raising parameter for the level-set function is explicitly computed to restore the target volume. The trick relies on an expansion of the volume enclosed by a level-set in terms of the height of this level-set:

|{φ -δc < 0}| := {φ<δc} 1dx = |{φ < 0}| + δc {φ=0} 1dσ + o(δc)
Therefore starting from a reference volume V 0 , one can define δc so that |{φ -δc < 0}| = V 0 + o(δc) by setting:

δc = V 0 -|{φ < 0}| |{φ = 0}| .
Regarding the conservation of the surface area, one approach, proposed in the context of vesicle or red blood cell simulations is to relax this conservation by introducing an area change energy with a high stiffness. In [9, 10] the authors showed that such an energy could be expressed thanks to ∇φ which records the area change of level-sets of φ when φ is advected by a divergence free vector field. By taking a high areal tension stiffness, ∇φ is kept close to 1 in a neighbourhood of the interface. In practice, this however induces a high stiffness of the numerical method, which leads to severe time step constraints. Another approach is to enforce exact zero surface divergence of the velocity field which advects the level-set function. This can be done using Lagrange multipliers [START_REF] Doyeux | Simulation of two-fluid flows using a finite element/level set method. application to bubbles and vesicle dynamics[END_REF][START_REF] Laadhari | Computing the dynamics of biomembranes by combining conservative level set and adaptive finite element methods[END_REF], leading to good area conservation, but at the expense of a bad conditioning of the underlying linear systems to be solve at each iterations, which increases the computational cost.

In this work, we introduce a correction in the spirit of Smolianski's trick to conserve both the volume and the area. From an optimisation point of view, our method can be regarded as a projection method, and the whole constrained Willmore flow as some kind of projected gradient method.

Raising a level-set function to achieve some given area and enclosed volume

Obviously, we cannot in general fulfill both the area and volume constraints by adding a constant to a level-set function. So let us look at the case where we would add a non-constant function δc : Ω → R to an arbitrary level-set function φ. Let H denote a one dimensional smoothed Heaviside function and ζ its derivative. A typical choice is:

H(r) =      1 2 r + 1 + 1 π sin(πr) |r| < 1 0 r ≤ -1 1 r ≥ 1 and ζ(r) = 1 2 (1 + cos(πr)) |r| < 1 0 |r| ≥ 1 (25)
We can then compute the enclosed volume of the raised level-set:

|{φ -δc < 0}| = lim ε→0 Ω H - φ -δc ε dx = lim ε→0 Ω H - φ ε dx + Ω 1 ε ζ φ ε δc dx + O δc 2 (26)
Note that if one adds to a level-set function a δc which averages to zero on {φ = 0}, the enclosed volume does not change.

Considering now the area, we have to take care that φ -δc is not in general a distance function:

|{φ -δc = 0}| = lim ε→0 Ω 1 ε ζ φ -δc ε |∇(φ -δc)|dx + O δc 2 = lim ε→0 Ω 1 ε ζ φ ε |∇φ|dx - Ω 1 ε 2 ζ φ ε |∇φ|δc dx - Ω 1 ε ζ φ ε ∇φ • ∇δc |∇φ| dx + O δc 2
Integrating by parts in the last term:

-

Ω 1 ε ζ φ ε ∇φ • ∇δc |∇φ| dx = Ω div 1 ε ζ φ ε ∇φ |∇φ| δc dx = Ω 1 ε 2 ζ φ ε ∇φ • ∇φ |∇φ| δc dx + Ω 1 ε ζ φ ε div ∇φ |∇φ| δc dx = Ω 1 ε 2 ζ φ ε |∇φ|δc dx + Ω 1 ε ζ φ ε κδc dx
where κ is the total curvature as introduced above. We finally get:

|{φ -δc = 0}| = |{φ = 0}| + lim ε→0 Ω 1 ε ζ φ ε κδc dx + O δc 2 (27)
Note that when φ is a distance function φ = d, the volume and area of the perturbed level-set d -δc simplify to formulae involving only "surface" integrals :

|{d -δc < 0}| = |{d < 0}| + {d=0} δc dσ + O δc 2
which boils down to the Smolianski formula when δc is constant, and

|{d -δc = 0}| = |{d = 0}| + {d=0} κδc dσ + O δc 2 .
Given a target volume V 0 and area A 0 , we propose the following fast method to correct a level-set function so that its zero level-set achieves those targets.

1. Raise φ by a constant -λ to achieve the right enclosed volume:

λ = V 0 -|{φ < 0}| Ω 1 ε ζ φ ε dx (28)
2. Raise φ = φ -λ by -µ(κ -κ), which does not change the enclosed volume, where the constant µ and κ are given by:

µ = A 0 -|{ φ = 0}| Ω 1 ε ζ φ ε κ(κ -κ) dx = A 0 -|{ φ = 0}| (κ 2 -κ 2 ) Ω 1 ε ζ φ ε dx κ = Ω 1 ε ζ φ ε κ Ω 1 ε ζ φ ε ( 29
)
Remark 3. In the practical numerical computations, the "thickness" ε of the interface is kept finite and proportional and close to the mesh size h: ε ∈ [h; 3h], following the usual level-set approach.

Clearly, the method fails if κ 2 = κ 2 , which is the case for a circle (resp. sphere). But in that case the perimeter (resp. area) and enclosed area (resp. volume) are linked and one cannot set one independently of the other.

Let us investigate how the former restoring of volume and area change the energy we are minimising. Assuming an energy of the form:

E[φ] = Ω Edν ε where dν ε = 1 ε ζ( φ ε
)dx and, with a gradient which can be written as

dE[φ](ψ) = Ω F ψdν ε . Proposition 1.
Let us consider a gradient-corrected method to minimise E iteratively by alternating steepest descent with corrections stages 1 and 2 described above. Then, at first order in the descent parameter, the energy does not increase during this iteration.

Proof. A steepest descent method to minimise this energy amounts to change φ to φ + ψ with ψ = -ρF , which decreases the energy, at first order, by ρ Ω F 2 dν ε . But this motion changed the enclosed volume and the area of interface. The volume is now |{φ + ψ < 0}| which correspond at first order to a change of ρ Ω F dν ε . Therefore the λ parameter computed above is equal to:

λ = ρ Ω F dν ε Ω dν ε =: ρ- Ω F dν ε
where from now on, we will denote with a dashed integral sign the mean value, i.e. the integral divided by the measure of the set (with respect to ν ε ) onto which the integral is taken:

-

Ω F dν ε = 1 ν ε (Ω) Ω F dν ε
The descent motion and this volume correction change the area by:

ρ Ω F κdν ε -λ Ω κdν ε = ρ Ω F (κ -κ)dν ε .
Thus the second step, area correction, computes a correction ν(κ -κ) with:

µ = ρ - Ω F (κ -κ)dν ε κ 2 -κ 2 .
Those two corrections raise the energy by:

λ Ω F dν ε + µ Ω F (H -H)dν ε = ρν ε (Ω) - Ω F dν ε 2 + - Ω F (H -H)dν ε 2 H 2 -H 2 .
Finally, the two successive corrections would not raise the energy more than it has been decreased by the descent method provided that the following inequality holds:

-

Ω F dν ε 2 + - Ω F (H -H)dν ε 2 H 2 -H 2 ≤ - Ω F 2 dν ε .
This is the object of the following lemma. We therefore justified that our descent-correction method does not increase energy.

Lemma 1. Let (X, Σ, µ) be a measure space, and A a measurable subset of X with µ(A) < +∞. Let f, g ∈ L 2 (A) with -A gdµ = 0. Then:

- A f dµ 2 - A g 2 dµ + - A f gdµ 2 ≤ - A f 2 dµ- A g 2 dµ ( 30
)
Proof. This is an easy extension of the Cauchy-Schwarz inequality. Indeed applying the latter for f --A f dµ and g we have:

- A f -- A f dx gdµ 2 ≤ - A f -- A f dµ 2 dµ- A g 2 dµ
Since g is of zero mean on A, we have for the left hand side of the former inequality,

- A f -- A f dµ gdµ = - A f gdµ
while concerning the right hand side we observe that:

- A f -- A f dµ 2 dµ + - A f dµ 2 = - A f 2 dµ
which leads to eq. ( 30).

Remark 4. In order to prevent undesired effects from outside of the neighbourhood of the interface we can slightly alter the method during the step recovering the area. We rescale with the function µ(κ -κ)

1 ε ζ φ ε rather than µ(κ -κ). This amounts to choose δc = (κ -κ) 1 ε ζ φ ε in eq. ( 27
) and to redefine κ (resp. H as:

κ = Ω 1 ε 2 ζ 2 φ ε κ Ω 1 ε 2 ζ 2 φ ε
In our simulations, we fix the value of the rescaling term µ(H -H) outside of the neighbourhood of the interface in order to have a neutral effect.

Remark 5. An alternative option is to rescale the volume and the area during the same step by solving the following system:

  Ω 1 ε ζ φ ε Ω 1 ε ζ φ ε H Ω 1 ε ζ φ ε H Ω 1 ε ζ φ ε H 2   • λ µ = V 0 -|{φ < 0}| A 0 -|{φ = 0}|
The first equation (resp. the second) corresponds to the volume expansion eq. ( 26) (the area expansion eq. ( 27)) with δc = λ + µH. It is important to note that the constant part µ H of the area term is absorbed inside λ. Both versions are equivalent in theory and yield similar numerical results.

Willmore flow with volume and area conservation constraints

Combining the algorithms alg. 3 and 5 for the Willmore flow derived in the previous section and the above method to recover volume and area, we can obtain an algorithm to compute the Willmore flow in both dimension 2 and 3 with conservation of volume and area.

The principle is to alternate one diffusion step moving according the Willmore flow and recovering the volume and area constraints before the redistanciation. Using proposition 1 we can ensure that the energy has globally decreased at the end of one whole step of the method.

The corresponding algorithms in 2D and 3D are described in alg. 6 and 7.

Numerical method and practical implementation

Finite Element library

The equations and integrals introduced above are solved or evaluated within a finite-element framework, using the Feel++-finite-element C++ library [START_REF] Prud | Feel++: A computational framework for galerkin methods and advanced numerical methods[END_REF]44], and in particular the LevelSet framework from the Feel++ toolboxes [42], which features a comprehensive and seamless parallel set of tools for this kind of surface capturing methods. This framework was preferred to a more classical finite difference approach for the simplicity of coding with high level templates that mimic the variational form of PDE. Moreover, our final goal is to solve fluid-structure interaction problems where vesicles are evolving inside a fluid which flows in a complex geometry like blood vessels.

More precisely, we use a continuous Galerkin variational approach, and discretise the resulting equations in space with Lagrange polynomials. Introducing T h ≡ {K e , 1 ≤ e ≤ N elt } a compatible tessellation of the computational domain Ω, and the corresponding discrete -unstructured -mesh Ω h = N elt e=1 K e , we define P k h ≡ P k h (Ω h ) as the finite-element space on Ω h spanned by Lagrange polynomials of order k.

Algorithm 6 2D Willmore flow with constant volume and area 1 while t < t f do Solve ∂ t φ -∆φ = 0 with initial condition

φ i = d (n) for times t = √ 2δt and t = δt/ √ → G √ 2δt , G δt/ √ 2 .
Compute the curvature H.

Compute D =        D E = -G √ 2δt + 2G δt/ √ 2 + H 3 2 δt 2 D CN = -G √ 2δt + 2G δt/ √ 2 + H 3 4 δt 2
Compute H(D) and its mean H(D).

Compute λ:

λ = V 0 -|{D < 0}| |{D = 0}| . Compute µ: µ = A 0 -|{D = 0}| Ω H(H -H) 1 ε ζ D ε dσ . Construct the new signed distance function d (n+1) = Redist D -λ -µ(H -H) . t ←    t + δt 2 (Euler) t + δt 2 2 (Crank-Nicolson)
Algorithm 7 3D Willmore flow with volume and area constraints 1 while t < t f do Solve ∂ t φ -∆φ = 0 with initial condition

φ i = d (n) for times t = √ 2δt and t = δt/ √ → G √ 2δt , G δt/ √ 2 . Solve ∂ t φ -∆φ = 0 with initial condition φ i = d (n)2 for times t = √ 2δt and t = δt/ √ → G d 2 √ 2δt , G d 2 δt/ √ 2 .
Compute the curvature H.

Compute D = Hd (n)2 2 + -G √ 2δt + 2G δt/ √ 2 • 1 -Hd (n) + -G d 2 √ 2δt + 2G d 2 δt/ √ 2 • H 2
Compute H(D) and its mean H(D).

Compute λ:

λ = V 0 -|{D < 0}| |{D = 0}| . Compute µ: µ = A 0 -|{D = 0}| Ω H(H -H) 1 ε ζ D ε dσ . Construct the new signed distance function d (n+1) = Redist D -λ -µ(H -H) . t ← t + 2δt 2 (Euler) t + δt 2 (Crank-Nicolson)

Solving the diffusion equation

From a numerical perspective, our Willmore flow alg. 3 and 5 mainly involve solving the diffusion equation ∂ t φ -∆φ = 0 with appropriate initial and boundary conditions. To this end, we use the P k h finite-element space discretisation introduced above, and a second order unconditionally stable Crank-Nicolson scheme for the time discretisation. The corresponding discrete variational problem then reads:

Find φ (n+1) ∈ P k h s.t. ∀ψ ∈ P k h , Ω h φ (n+1) δt ψ + 1 2 ∇φ (n+1) • ∇ψ = Ω h φ (n) δt ψ - 1 2 ∇φ (n) • ∇ψ + ∂Ω h N • ∇φ (n) ψ (31)
where the superscript indices denote the time iterations. Note that the last term in eq. ( 31) comes from the integration by parts of the diffusive terms, and is somehow similar to an explicit discrete Neumann boundary condition preserving the normal gradient of φ, namely N •∇φ (n+1) = N •∇φ (n) with N the exterior normal at the boundary of Ω h . In practice, the best choice of boundary conditions for our distance diffusion problem is an open question, since no natural condition emerges from our analysis, which focuses on the neighborhood of the interface Γ, far from the domain boundary. As pointed out in remark 1, we would also like to use in some situations an Euler scheme, which is low order and diffusive, which may be more adapted in this context than a dispersive scheme. In that case one should directly use the discrete expansion, which differs at order two (since the scheme is order one):

φ n+1 -δt∆φ n+1 = φ n indeed gives φ n+1 = (id -δt∆) -1 (φ n ) = φ n + δt∆φ n + δt 2 ∆ 2 φ n + o(δt 2 ) (32)
which differs from of ( 5) by a factor 2 in the second order term. The corresponding discrete variational problem then simplifies to:

Find φ (n+1) ∈ P k h s.t. ∀ψ ∈ P k h , Ω h φ (n+1) δt ψ + ∇φ (n+1) • ∇ψ = Ω h φ (n) δt ψ + ∂Ω h N • ∇φ (n) ψ (33)
In our numerical results, both schemes give qualitatively the same equilibrium shapes. However the Crank-Nicolson scheme seems more accurate, while the Euler scheme was preferred in dimension 3 where its diffusive behavior brings more stability and smoothness of the interface.

Computation of the curvature

Our algorithms for the Willmore flow in both 2D alg. 3 and 3D alg. 5 require the computation of the curvature of the surface Γ in the correction term. We first present the direct methods currently used in Feel++. We then introduce our diffusion methods to compute the curvature, which are essentially the classical convolution/thresholding algorithms of order 1 and 2. Using the direct methods as references, we discuss the efficiency of our method and focus on the choice of the time step selection for the diffusion steps.

Direct methods

From a usual level-set perspective, the curvature of Γ can be computed directly using the divergence of the level-set function ∆d. For low-order discretisations however, this requires the use of specific strategies, as standard finite-element derivation decreases the discretisation polynomial order by 1, which forbids derivations of degrees higher than the polynomial order. To circumvent this issue, we use a classical Galerkin projection in L 2 (Ω h ) to maintain the element in P 1 h . Namely the computation of the gradient ∇d is as following:

Find g ∈ P 1 h s.t. ∀v ∈ P 1 h , Ω h g • v = Ω h ∇d • v
We can then obtain the curvature with two successive derivation/projection steps by computing the divergence of the gradient of the distance function:

(1) :

Find g ∈ P 1 h s.t. ∀v ∈ P 1 h , Ω h g • v = Ω h ∇d • v (2) : Find H ∈ P 1 h s.t. ∀w ∈ P 1 h , Ω h Hw = Ω h div g w (34)
However the L 2 projection fills the missing information of the lower order elements with noisy values, so that one usually resort to a smoothed L 2 projection method, where a small diffusion term -η∆H is added:

(1) :

Find g ∈ P 1 h s.t. ∀v ∈ P 1 h , Ω h g • v = Ω h ∇d • v (2) : Find H ∈ P 1 h s.t. ∀w ∈ P 1 h , Ω h Hw + η Ω h ∇H • ∇w -η ∂Ω h ∇H • N w = Ω h div g w. ( 35
)
where the smoothing coefficient η is typically taken equals to 0.03h. We will refer to eq. ( 34) as the L 2 projection and to eq. ( 35) as the Smoothed projection method.

Using the diffusion of the signed distance function

As mentioned above, we can also retrieve the curvature of the surface represented with the level-set using the diffusion of the signed distance function. Recalling the analytical expansion of the diffusion solution eq. ( 6), we observe that the curvature of the surface shows up in the first order term in δt, and can thus be retrieved as:

H = G δt -d (n -1)δt + O(δt). ( 36
)
As explained in section 2.3, at the expense of an additional resolution of the diffusion equation, we can improve the order of accuracy and compute the combination eq. ( 16) which eliminates the second order term. We can then compute H at order δt 2 as:

H = a 1 G α + b 1 G β -3 √ 2 2 d (n -1)δt + O(δt 2 ) (37)
In the following we shall refer to the O(δt) method eq. ( 36) as the Order 1 diffusion and the O(δt 2 ) method eq. (37) as the Order 2 diffusion. Keep in mind however that these approximation orders refer to the auxiliary diffusion time-step and not to the usual space discretisation order.

Before comparing these diffusion methods to the direct ones mentioned in section 4.3.1, we shall first study in more detail the role of the diffusion time-step. Figure 1 shows the evolution of the L 2 error of the curvature computed with both diffusion methods as a function of the diffusion time-step. The numerical test was performed on a two-dimensional circular surface of curvature 1 with a mesh size h ≈ 0.02, and the error was computed as

err(H; Γ) = Ω (H -H th ) 2 1 ε ζ φ ε Ω 1 ε ζ φ ε , (38)
where H th stands for the theoretical value of H. In the case of a circle, it amounts to the inverse of the radius. The spreading parameter ε was chosen for all the simulations so that the interface span about 2 mesh cells, i.e. ε ≈ h. 36) and Order 2 eq. ( 37) diffusion methods as a function of the diffusion time-step δt. We observe the expected convergence orders as the time-step decreases, and a saturation of the error for low time-steps which is related to space discretisation error terms.

As expected, we observe convergence of both methods at order O(δt) and O(δt 2 ) respectively, as well as a saturation of the error for small time-steps, which is due to the spatial discretisation errors. We also note that the Order 2 diffusion method features an optimal error for larger time-steps, which can prove useful when using the same time-step for the flow and the diffusions (c.f. remark 2).

These numerical observations, as well as the analytic expansion of G δt also suggest to choose the time-step following the heuristic strategy:

Hδt ≈ h (39)
with H an a priori estimate of the typical curvature of the surface. From a numerical point of view, this corresponds to having the diffusion move the d = 0 level-set a few mesh elements. This heuristic can also be adapted to the case of geometrical flow simulations, and naturally formulates an adaptive time-stepping method, as presented in section 4.5.

Comparison of the different methods

We now turn to the comparison of the "direct" and "diffusion" curvature methods presented above. We again use a fixed two-dimensional circular interface of radius 1, and vary the characteristic mesh size h used to compute the L 2 error on the curvature estimate along the surface as a function of h for the four methods. We set the smoothing parameters η = 0.03h and the diffusion time-steps as δt = 0.7h. The resulting errors are plotted in log-log scale in fig. 2.

We find that the Order 1 and Order 2 diffusion methods behave similarly, and compare accurately with the Smoothed projection method. These three methods give O(h) convergence, which is satisfactory for our low-order discretisation. As anticipated, the L 2 projection method does not seem to converge.

In our simulations, we shall therefore use the Order 2 diffusion method, as the two diffusion equation resolutions for G α and G β are needed anyway for the computation of the Willmore flow, and this method provides the best curvature estimate at no additional cost.

Level-set redistanciation

The redistanciation step is a crucial aspect of our algorithm, as it allows the iteration to proceed while preserving the correct expansions of the diffusion solutions close to the interface. This step can moreover become the main bottleneck of the whole algorithm in terms of both efficiency and stability if not processed carefully. From a general point of view, performing the Redist step in alg. 3 and 5 amounts to computing the distance to an interface Γ located by the 0-level of some function. Stated differently, we need to solve the eikonal equation:

|∇φ| = 1, φ(Γ) = 0 (40)
which is a challenging boundary-value non-linear equation which has received much attention since the work of Sethian [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF] who proposed an iterative upwind fast-marching algorithm to solve the eikonal equation starting from the interface and propagating the information outward.

In our finite-element framework, we use a parallel fast-marching algorithm inspired from [62] but adapted to arbitrary -possibly unstructured -meshes. The local -element-wise -eikonal equations are solved with a QR decomposition, and the fast-marching alternates between local-domain solves and ghost inter-domain updates until global convergence is obtained.

In order to ensure good stability properties of our method, and in particular to prevent the fast-marching method from introducing spurious motion of the interface, the initialisation of the algorithm near the implicitly defined interface is crucial. To this end, we rescale the level-set function φ by 1 |∇φ| on the elements which intersect the interface, which exactly solves the local eikonal equation when φ is a piecewise linear polynomial (φ ∈ P 1 h ), before applying the fast-marching algorithm both inward (φ < 0) and outward (φ > 0) starting from the values encompassing the interface.

Adaptive time-stepping

As can be seen from eq. ( 17), given some δt in the resolution of the heat equations, our diffusionredistanciation schemes actually solves the corresponding Willmore flow with time-step δt 2 , namely

d (n+1) = d (n) -W δt 2 .
The flow time-step is therefore not constrained by any stability issue, but controls the accuracy of the solution, since the expansions performed in section 2.1 are valid up to terms of order O(δt 3 ).

The O(h) errors pertaining to the computation of φ (n+1) in eq. (31) using P 1 h Lagrange elements and the expansion (6) however suggest to use a time-step like δt ∼ h κ in the heat equation. This obviously requires the a priori computation of the interface curvature, which is not desirable in our framework. Instead, we therefore use an explicit adaptive time-stepping strategy by setting

δt (n+1) = min Γ h κ (n) (41)
where κ (n) is the curvature at the previous iteration and h is a measure of the mesh size, here considered for elements crossed by the interface -e.g.the minimal diameter of the elements crossed by the {φ (n) = 0} level-set.

Numerical illustrations

In this section, we illustrate the efficiency our methods by applying them to two classical applications. First, we apply our 3D algorithm for the Willmore flow alg. 5 to the Willmore problem. Then we apply our algorithm for Willmore flow with area and volume conservation alg. 7 to the computation of an equilibrium shape of a red blood cell.

Willmore problem

The Willmore problem refers to the study of the minimisers of the Willmore energy function E W . The case of compact surfaces of genus 0 is trivial as the minimisers are the spheres -recall that the Willmore energy being scale invariant, all the spheres have the same energy.

In the case of compact surfaces of genus 1, the Clifford torus, the torus with a ratio √ 2 between its radii, minimises the Willmore energy. The proof of this conjecture made by Willmore [60] in 1965 has been established recently in [39]. The conjectures for higher genus order [START_REF] Hsu | Minimizing the squared mean curvature integral for surfaces in space forms[END_REF] are still to be proven.

In the following, we present 2D and 3D simulations of such simple surfaces which evolve according to the Willmore flow using our diffusion-redistanciation algorithms. The existence of analytic minimisers will then provide solid means to evaluate our numerical approach.

2D Willmore flow of a circle

To assess the accuracy of our diffusion-redistanciation scheme, we first simulate an initial two-dimensional circular interface evolving according to the Willmore flow, which can also be computed analytically. As it evolves with the flow, the interface should remains circular, and increase it radius r as r(t) = (r 0 + 2 t) 1/4 (42) with r 0 the initial radius. Note that this law, which comes from the analytic flow equation dr dt = H 3 2 = 12r 3 provides a good test of our algorithm, since the right hand side is actually obtained through a delicate compensation between the second-order term w (c.f. ( 10)) in the diffusion of d and the correction term c obtained with the diffusion of d 2 (c.f. ( 12)).

Our simulation was run with r 0 = 1 in a square domain of side length 6 with an unstructured mesh of typical size h ≈ 0.02 and the adaptive time-step strategy presented in section 4. [START_REF] Biben | Tumbling of vesicles under shear flow within an advected-field approach[END_REF].

Figure 3 shows some snapshots of the simulation and the evolution of the circle radius. It illustrates the numerical accuracy of our scheme, which succeeds in preserving the symmetry of the shape as it evolves, and compares quantitatively with the exact solution. Figure 4 shows the error on the computed radius of the circle as compared to the theoretical value (42) at t = 3 for different mesh sizes. The dependence of this error on the mesh size demonstrates the convergence of our method. Note however that the fitted convergence order -≈ 1.72 -is not strictly equivalent to the usual finite-element convergence order, as it incorporates effects related to the dependence of the time step on the mesh size, which is chosen following eq. ( 41) as the role played by the time step in the derivation of our algorithm naturally imposes such spatio-temporal coupling. The linear dependence of the time-step with respect to the mesh size would suggest that the actual spatial convergence order is ∼ 0.7 -0.8, which also seems supported by the energy convergence shown in fig. 20. Figure 4: Evolution of the error err(r(t = 3)) ≡ |r h (t = 3)r(t = 3)| of the radius of the circle at t = 3 as a function of the mesh size. The figure also shows the best mean-square fit of the data by a power function 14.48 h 1.72 . Note however that the displayed slope is not strictly equivalent to the usual convergence order, as discussed above.

3D Willmore flow of a torus

The convergence of a torus towards a Clifford torus, whose ratio between major and minor radii equals √ 2, is a good test for our numerical approach, as it challenges both the accuracy of the numerical flow and the stability of the algorithm as it approaches the minimising surface, which features a rather small hole, which can easily be filled by an inaccurate flow to minimise the energy further, as the Willmore energy of any sphere -4 π -falls below the minimal one for tori, which is 2 π 2 .

Simulations have been performed using a phase-field model in [6] but the results are only qualitative. To the extend of our knowledge, there is no other published numerical work to compare with our results.

Figures 5 and6 show the results of our three-dimensional Willmore flow algorithm 5 starting from a torus with major a = 2 and minor b = 0.5 radii. The simulation was run in a cuboid with lengths 6 × 4 × 6 and an unstructured mesh of typical size h ≈ 0.02. We observe that the surface flows to the expected Clifford torus, and that the estimated resulting energy and radii ratio are in good agreement with the expected values for such a torus. However we were unable to stabilise the optimal shape once it reaches the Clifford torus. We observe in Figure 5 the energy fluctuating and even increasing slowly beyond t = 3, while in Figure 6 the ratio a/b does not seem to converge. The reason behind this behaviour is that we are computing a and b from the surface and enclosed volume assuming that we are in the class of toroidal surfaces. However it turns out that our algorithm, when reaching the Clifford torus, tries to minimise the energy by slipping out of this class (we remind that we are making a full 3D computation not imposing any symmetry). There the central and peripheral circles have no longer exactly the same center. It seems that while reaching the Clifford torus, the energy is flattening and it is more likely that discretisation errors drive the surface out of the tori class. Once this symmetry breaking has occurred, the hole is then quickly closed and we fall in the 0 genus surfaces class of lower energies. This raises the question on how to constraint somehow the surface to stay in the tori class (without doing axisymmetric computations). For now we did not find any simple algorithm to tackle that problem, which while interesting, is not directly related to our aimed application to red blood cells.

Vesicles equilibrium shapes

We now consider the simulation of vesicles, which provide the most popular and simple model for capsules, closed thin shell or cells, like human red blood cells. From a physical point of view, the model consists in an elastic thin shell of fixed area enclosing some incompressible fluid, so that the inner volume is also fixed. As such, the vesicle only deforms through bending of the shell, and is thus controlled by a Willmore-like -or Canham-Helfrich energy [START_REF] Helfrich | Elastic properties of lipid bilayers: theory and possible experiments[END_REF][START_REF] Peter B Canham | The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell[END_REF] (see e.g. [51] for a presentation of the model in a biological context), which reduces to the standard Willmore energy when the membrane has no spontaneous curvature -the unconstrained (open) membrane is flat, which we assume in this work.

In the following, we shall use our diffusion-redistanciation algorithms with the constraints of constant volume and area of the evolving surface to compute the 2D and 3D equilibrium shapes of vesicles with various reduced volumes. We show in the following the results of the "Euler" variants of alg. 6 and alg. 7. Figure 6: Willmore energy and geometrical quantities of the torus. We observe that the flow decreases the energy to about 20.8, which is within 5% from the theoretical Clifford torus energy 2π 2 ≈ 19.74. We also observe that at equilibrium, the ratio between the radii is indeed about √ 2. Beyond t = 3, the surface leaves the toroidal class, and the computed a and b are not anymore relevant (see above).

We have also performed the same simulations using the "Crank-Nicolson" variants of the algorithms, which yielded very similar shapes and results, but featured some small -stable -oscillations at the equilibrium, as could be expected from a non-diffusive discretisation scheme. The resulting equilibrium shapes are compared to the ones from the literature when available, and assessed from energetical and stability points of view to illustrate the good properties of our numerical method.

Computation of the equilibrium shapes of 2D vesicles

We first consider the two-dimensional case, and compute the Willmore flow under the constraint of constant inner volume and surface area of an initial ellipse with semi-minor and -major axes adjusted according to a chosen reduced volume

ν = 4πV A 2 = V π A 2π 2 (43)
The reduced volume represents the ratio between the volume of the ellipse and the volume of a circle of same area. Note that the conservation of both the volume and area of the vesicle naturally entails that the reduced volume is also preserved. The simulations were run in a square domain [-4, 4] 2 with a structured triangle mesh with mesh size h = 0.04 and adaptive time-step strategy. We illustrate the convergence of our method to an equilibrium shape in fig. 7 and compare our results with the ones obtained using a lattice-Boltzmann method in [26] in fig. 8. We observe that the shapes computed with our algorithm are in excellent agreement with the ones obtained from direct lattice-Boltzmann simulations. We also illustrate the good volume and area conservation properties of our algorithm in fig. 9. The evolution of the Willmore energy along the flow is shown in fig. 10.

Computation of vesicles equilibrium shapes in 3D

The landscape of vesicle equilibrium shapes is much more rich in three dimensions, as the Gauss curvature comes into play and allows for multiple local minima for shapes with the same topological class and invariants. Axisymmetric equilibrium shapes have been mapped out in a phase diagram in [START_REF] Seifert | Shape transformations of vesicles: Phase diagram for spontaneouscurvature and bilayer-coupling models[END_REF], and some non-axisymmetric shapes have been computed in [25, [START_REF] Du | Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions[END_REF][START_REF] Du | Adaptive finite element method for a phase field bending elasticity model of vesicle membrane deformations[END_REF]6], but the general phase diagram for Willmore energy minimisers is still an open question.

In our case, we consider the evolution of ellipsoids under the constant volume and area Willmore flow. As such, we expect to recover the prolate and oblate axisymmetric shapes of zero genus and spontaneous curvature referenced in [START_REF] Seifert | Shape transformations of vesicles: Phase diagram for spontaneouscurvature and bilayer-coupling models[END_REF] depending on the initial -conserved -reduced volume,

ν = 6 √ πV A 3 2 = V 4 3 π A 4π 3 2 . ( 44
)
As shown in [START_REF] Seifert | Shape transformations of vesicles: Phase diagram for spontaneouscurvature and bilayer-coupling models[END_REF], the prolate and oblate shapes are both local minimisers for ν ≈ 0.51. While the prolate [26] using a lattice-Boltzmann method for reduced volumes ν = 0.9, 0.8, 0.7 and 0.6. We also give our result for a reduced volume ν = 0.5. shapes exist for all possible ν, the oblate ones self-intersect below ν ≈ 0.51, and the stomatocyte shapes become the only non-prolate feasible minimisers. In this work, we restrict ourselves to ν ≥ 0.6, and study both the oblate and prolate cases, as our algorithm appears able to capture the local minimising shapes in a stable way. Fewer quantitative works are available for comparison in this three dimension case. We can mention [4, [START_REF] Du | Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions[END_REF]20,[START_REF] Du | Adaptive finite element method for a phase field bending elasticity model of vesicle membrane deformations[END_REF] where a mix of qualitative and quantitative results are presented. In [4, 14], a phase field method was developed and qualitative results were given. In [START_REF] Du | Adaptive finite element method for a phase field bending elasticity model of vesicle membrane deformations[END_REF] an adaptive version of the phase field model was proposed with some quantitative results in terms of energy, which are hard to compare with theoretical known values of Seifert [START_REF] Seifert | Shape transformations of vesicles: Phase diagram for spontaneouscurvature and bilayer-coupling models[END_REF]. The approach of Feng and Klug [20] is based on surface finite element and is more quantitative in its results as far as the energy is concerned.

Oblate case. The oblate case simulations were run in a cuboid with lengths 3.6 × 6.4 × 6.4 and mesh size h ≈ 0.04. The initial level-set functions were taken as signed distance functions to oblate ellipsoids with semi-minor axis a and semi-major axes b = c adjusted to get a volume 4π 3 for all the simulations and a corresponding reduced volume ν = 0.6, 0.65, 0.7, 0.8 or 0.9.

Figure 11 shows the equilibrium shapes obtain with our constrained diffusion-redistanciation algorithm 7 for the different reduced volumes. We also show two-dimensional cuts of these equilibrium shapes in fig. 12 to illustrate the good symmetry conservation properties of our algorithm. Figures 13 and14 show the evolution of the Willmore energy, and the area and enclosed volume of the surface as it evolves. Prolate case. The prolate case simulations were run in a cuboid with lengths 9×4×4 and mesh size h ≈ 0.04. The initial level-set functions were chosen as in the oblate case, but using prolate ellipsoids. The resulting equilibrium shapes and cuts are shown in figs. 15 and 16 respectively, while figs. 17 and 18 show the evolution of the surface energy, area and enclosed volume.

. We can observe for all the -oblate and prolate -three-dimensional vesicle simulations that the computed equilibrium shapes are axisymmetric as expected, even though our algorithm does not enforce this property. We also note that our algorithm seems very robust in this case from the energy minimisation point of view, as it reaches the minimum rather quickly and then preserves it for a long time without any numerical artifact.

To compare our results quantitatively, we plot in fig. 19 the equilibrium energies of our vesicles together with the ones obtained in [START_REF] Seifert | Shape transformations of vesicles: Phase diagram for spontaneouscurvature and bilayer-coupling models[END_REF] and [20]. We also highlight the numerical convergence of our method in fig. 20 where we plot the relative error of the Willmore energy of our final equilibrium shape (for the reduced volume ν = 0.65) as a function of the mesh size. Our results seem in good agreement with the ones obtained with direct energy minimisation or meshed surface evolution, despite a small ∼ 10% overestimate for the oblate smallest reduced volume cases, which seem related to typical finite-element numerical errors. As shown in fig. 20, the equilibrium shape energy converges to the expected value (taken from [START_REF] Seifert | Shape transformations of vesicles: Phase diagram for spontaneouscurvature and bilayer-coupling models[END_REF] for this axisymmetric case) as h 0.84 , which is rather positive for low-order (P 1 ) simulations of such high-order effects.

Conclusion

We propose in this article a numerical method to predict the position of an interface in 2D or 3D moving according to the gradient of an higher order energy such as the squared mean curvature (Willmore flow) and subject to constraint of area and enclosed volume conservation. This problem is of interest to devise semi-implicit schemes for fluid-structure solvers where the full interaction between immersed vesicles and fluids are involved. As a first step, the present work was restricted to the motion of vesicles membranes minimising their mean curvature at fixed surface area and enclosed volume. The test case considered corresponds to reach an equilibrium shape starting from an ellipsoid, depending on its closeness to a sphere or more elongated shape. Even in this delimited setting, the problem is tricky since the constraints are highly nonlinear and the energy involved of higher order. Classical approaches through level-set or phase field methods lead to fourth order PDE to solve, while purely Lagrangian methods dealing with a surface mesh are not well suited to be included in a fluid-structure coupling procedure due to the interpolations required between Eulerian and Lagrangian representations. In this work, we build a method where only heat equations are solved to compute the right flow of the mean curvature energy.

Based on diffusion-redistanciation schemes introduced in [START_REF] Esedoglu | Diffusion generated motion using signed distance functions[END_REF], our first contribution was to provide a more intrinsic formulation of these methods, which opens the way to study more easily other kind of higher order energy. While this study focused on the mean curvature flows, one could for instance also consider problems where the Gaussian curvature is involved. We also extended our geometrical flows algorithms to the evolution of surfaces with conserved area and enclosed volume. The originality of our approach relies on a formula providing explicitly the projection of the unconstrained Willmore flow on motion conserving area and enclosed volume. This is a big advantage in comparison to other methods where the area constraint is penalised (and therefore not exactly fulfilled or leading to stiff problems) or non-linearly enforced (at a high computational cost). W /E th W with E th W taken from [START_REF] Seifert | Shape transformations of vesicles: Phase diagram for spontaneouscurvature and bilayer-coupling models[END_REF], and the mesh sizes range from h ≈ 0.02 to h ≈ 0.06. The plot is in log-log scale.

Our first test case was devoted to the torus, where we observed convergence to the Clifford torus as expected. Our methodology enjoys numerical convergence and compares well with existing numerical results to compute equilibrium shapes of vesicles in dimension 2. In dimension 3, while our results are in good agreement for prolate forms, some discrepancy occurs for oblate vesicle shapes. This is due to the low order of the finite element method used (P 1 element), in contrast with the C 1 surface finite element representation of [20]. However we observe numerical convergence in dimension 3 as well. The use of polynomial discretisation with higher regularity (such as Hermite elements) would likely cure that problem without any further modification of our algorithm. One remarkable feature of our algorithm is to provide consistent results without any mesh refinement strategy. This is of paramount importance for its possible use in a full threedimensional fluid-structure problems, where such remeshing could very quickly lead, for realistic situations, to intractable computational complexity.

Another big advantage of our approach is to be very easy to implement: high order geometric quantities are computed from diffusion equation of a distance function, which is easily handled by any finite element library.

In a forthcoming work, we will use such a predictor of motion to devise a semi-implicit scheme for the full fluid-membrane coupling problem.

Figure 1 :

 1 Figure1: L 2 error on curvature estimates using the Order 1 eq. (36) and Order 2 eq. (37) diffusion methods as a function of the diffusion time-step δt. We observe the expected convergence orders as the time-step decreases, and a saturation of the error for low time-steps which is related to space discretisation error terms.

Figure 2 :

 2 Figure2: Evolution of the L 2 error as a function of the mesh size for the different methods. The L 2 projection shows no evident convergence while the Order 1 diffusion, the Order 2 diffusion and the Smoothed projection methods yield linear convergence.

  Snapshots of the interface at times t = 0, 1, 2 and 3. Evolution of the radius of the circle: comparison with the exact solution.

Figure 3 :

 3 Figure 3: 2D Willmore flow of a circle.

Figure 5 :

 5 Figure 5: 3D Willmore flow of a torus toward the Clifford torus. The evolving surface is colored with the value of the curvature while the Clifford torus is represented in semi-transparent white.

  Evolution of the Willmore energy E W = Γ H 2 . Evolution of the torus major and minor radii.

Figure 7 :

 7 Figure7: Convergence of alg. 6 to an equilibrium shape for the ν = 0.5 case, starting from an ellipse with semi-major and -minor axes 2.09021 and 0.478421. The curves are colored with the value of the curvature.

5 Figure 8 :

 58 Figure8: Comparison of the 2D equilibrium shapes obtained by alg. 6 with the ones obtained by[26] using a lattice-Boltzmann method for reduced volumes ν = 0.9, 0.8, 0.7 and 0.6. We also give our result for a reduced volume ν = 0.5.

5 Figure 9 : 5 Figure 10 :

 59510 Figure9: Evolution of the 2D vesicle volume and "area" (perimeter in this case) with the volume-and area-preserving alg. 6. The algorithm displays very good conservation properties : the relative volume and area changes are less than respectively 1 • 10 -4 and 5 • 10 -4 for the most difficult ν = 0.5 case, and the final equilibrium shapes volumes and areas relative changes are 1 • 10 -5 for all the cases.

Figure 11 :

 11 Figure11: 3D oblate vesicles equilibrium shapes obtained with alg. 7 for the reduced volumes ν = 0.9, 0.8, 0.7, 0.65 and 0.6.

Figure 13 :

 13 Figure13: Evolution of the 3D oblate vesicle volume and area with the volume-and area-preserving alg. 7. The 3D algorithm also displays good conservation properties: the relative volume and area changes are respectively 5 • 10 -5 and 2 • 10 -4 for all the cases.

Figure 14 :

 14 Figure 14: Evolution of the Willmore energy E W = Γ H 2 (18) of the 3D oblate vesicle with the volume-and area-preserving alg. 7.

Figure 15 :x -y slice x -z slice y -z slice Figure 16 :

 15slice16 Figure 15: 3D prolate vesicles equilibrium shapes obtained with alg. 7 for the reduced volumes ν = 0.9, 0.8, 0.7, 0.65 and 0.6.

Figure 17 :Figure 18 :

 1718 Figure 17: Evolution of the 3D prolate vesicle volume and area with the volume-and area-preserving alg. 7. The 3D algorithm also displays good conservation properties: the relative volume and area changes are respectively 5 • 10 -6 and 3 • 10 -5 for all the cases.

Figure 19 :Figure 20 :

 1920 Figure19: Willmore energy of the 3D vesicle equilibrium shapes as a function of the reduced volume, in the oblate and prolate cases. We compare our results to axisymmetric energy minimisation (shooting method) from[START_REF] Seifert | Shape transformations of vesicles: Phase diagram for spontaneouscurvature and bilayer-coupling models[END_REF] and to the "subdivision thin shell" vesicle simulations from[20].

} .

Acknowledgements

Most of the computations presented in this paper were performed using the Froggy platform of the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche, and the Atlas cluster from the Research Institute in Advanced Mathematics (IRMA -UMR7501).