Computing arbitrary Lagrangian Eulerian maps
for evolving surfaces
Abstract
The good mesh quality of a discretized closed evolving surface is often compromised during time evolution. In recent years this phenomenon has been theoretically addressed in a few ways, one of them uses arbitrary Lagrangian Eulerian (ALE) maps. However, the numerical computation of such maps still remained an unsolved problem in the literature. An approach, using differential algebraic problems, is proposed here to numerically compute an arbitrary Lagrangian Eulerian map, which preserves the mesh properties over time. The ALE velocity is obtained from a simple spring system, based on the connectivity of the nodes in the mesh. We also consider the algorithmic question of constructing acute surface meshes. We present various numerical experiments illustrating the good properties of the obtained meshes and the low computational cost of the proposed approach.
Keywords: arbitrary Lagrangian Eulerian map, evolving surfaces, evolving surface PDE, ESFEM
AMS: 65M50, 34A09, 35R01
1 Introduction
Partial differential equations on evolving surfaces with a given velocity have been discretized using a huge variety of methods. Probably one of the most popular is the evolving surface finite element method developed by Dziuk and Elliott in [7].
As it was pointed out by Dziuk and Elliott early on, in Section 7.2 of [7]: ”A drawback of our method is the possibility of degenerating grids. The prescribed velocity may lead to the effect, that the triangulation is distorted”. The same issue occurs for problems in moving domains.
In recent years this problem has been theoretically addressed in a few ways:

The ALE version of the evolving surface finite element method has been proposed by Elliott and Styles in [13], where a better triangulation is obtained using an ALE map, i.e., allowing the nodes of the mesh to evolve with a velocity having an additional tangential component compared to the (pure normal) surface velocity.
We propose here to compute an arbitrary Lagrangian Eulerian map for closed evolving surfaces by integrating a differential algebraic equation (DAE) system for the nodes. We use a not necessarily tangential ALE velocity to achieve good mesh quality, while enforcing the points to stay on the surface.
To our knowledge there is no such ALE algorithm available in the literature, in contrast with the many papers on the theory of numerical methods involving ALE maps for both closed evolving surfaces and moving domains.
Many experiments have been presented in the above references, especially see [13, 14, 25], where smaller discretization errors have been obtained by solving evolving surface problems on ALE meshes. The ALE maps used in these experiments were unrealistic, obtained analytically from an a priori knowledge on the surface and its evolution, using deep understanding and structure of the signed distance function (which defines the surface). No general ideas on the computation of ALE maps for evolving surfaces have been proposed in these papers. Numerical analysis of the ALE evolving surface finite element method has been studied in [14] and [25].
For nonALE methods it is usual that nodes disappear and/or new nodes are added to the mesh (see, e.g., [11, Figure 14]). From a theoretical point of view this seems a minor issue, however from the implementation side this is undesirable, since then a map – between the old and the new mesh – has to be computed after each remeshing process. In most cases this is a nontrivial and costly task, and hence, should be avoided.
In the present paper we propose a general algorithm to compute a suitable ALE map, without any a priori knowledge, for meshes of closed evolving surfaces.
The approach is based on the following idea: Usually, the surface evolution is given by an ODE system with a surface velocity assumed to be normal. We use an additional tangential velocity for a possibly degenerated mesh to improve grid quality. In general such tangential velocities are not straightforward to construct, therefore we use a not necessarily tangential ALE velocity and introduce a constraint to keep the nodes on the surface. Altogether this is finally formulated as a DAE system for the nodes. The ALE velocity is based on a spring system, where the nodes are connected along the edges by springs. The numerical solution of the DAE system gives the new mesh. We use implicit Runge–Kutta methods (in particular Radau IIA methods), and a more efficient splitting method, combined with explicit Runge–Kutta methods, to integrate the system in time.
The computation of the arbitrary Lagrangian Eulerian mesh here is free of any a priori knowledge in the following sense: the algorithm uses the distance function at each time, but it does not use its structure or any other special properties of it, unlike the ALE maps from the literature.
This approach for closed evolving surfaces can be used as a tool in the computation of ALE meshes for moving domains: In [17, Section 2.4] arbitrary Lagrangian Eulerian maps for moving domains are obtained by solving a parabolic problem, or the corresponding stationary problem, while in [15] an elastodynamic equation system is used for the same purpose. However, for these approaches the evolution of the boundary still needs to be known a priori. The problem of numerically finding such a boundary evolution has not been solved in these papers. In fact, to determine such a boundary evolution is equivalent to finding an ALE map for a lower dimensional closed surface, which is the same problem as we consider in this paper. Hence, the algorithm proposed here can also serve as a tool to compute boundary evolutions, which can be used together with the well understood classical ALE methods for moving domains, for instance the ones proposed in [15, 17].
We give some further details on possible extensions of the proposed algorithm: to handle other mesh properties (e.g., acuteness), an adaptive version, and a local version as well.
We present various numerical experiments illustrating the validity of the differential algebraic model, and also the performance of the proposed algorithm compared to the ALE maps given in the literature. We also report on computational times in the case of evolving surface PDEs.
2 Evolving surfaces, ALE maps and PDEs
As our main motivation lies in the numerical solution of parabolic PDEs on evolving surfaces we shortly recap the setting of [7]. We will also use this setting as an illustrative background to the proposed algorithm.
Let , , be a smooth evolving closed hypersurface. Further, let the evolution of the surface given by the smooth velocity , assumed to be normal. Let denote the material derivative of , the tangential gradient is denoted by and given by , with unit outward normal . We denote by the Laplace–Beltrami operator.
We consider the following linear evolving surface PDE:
on  (1)  
on 
Basic and detailed references on evolving surface PDEs and on the evolving surface finite element method (ESFEM) are [7, 8, 9] and [27]. The ALE version of the evolving surface finite element method has been proposed in the paper [13], which also contains a detailed description and many experiments. While numerical analysis of full discretizations can be found in [14] and [25], the former is more concerned about spatial discretizations and BDF methods of order 1 and 2, while the latter is more focused on highorder BDF and Runge–Kutta time discretizations.
As we aim at the numerical solution of the PDE (1) using the surface finite elements developed by Dziuk and Elliott [7], the surface is discretized using a triangular mesh. The description of representation of evolving surfaces and that of discrete surfaces can be found in the following subsections.
2.1 Surface representations
Evolving surfaces are usually described in two ways, which have different advantages, hence we will use both of them for various purposes.
Distance function representation. Based on a signed distance function the evolving dimensional closed surface is given by
with a function (whose regularity depends on the smoothness of the surface), cf. [5, 7].
Diffeomorphic parametrization. The surface can also be described by a diffeomorphic parametrization, cf. [24] and [7].
We consider the evolving dimensional closed surface as the image
of a sufficiently regular vectorvalued function , where is the smooth closed initial surface, and . 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 parametrisation also satisfies the ODE system, for a point ,
(2) 
where is the velocity of the surface. Note that for a known velocity field , the position at time of the particle with label is usually obtained by solving the ordinary differential equation (2) from to for a fixed .
2.2 Surface approximation
The smooth initial surface is approximated by a triangulated surface , i.e., an admissible family of triangulations of maximal element diameter ; see [7] for the notion of an admissible triangulation, which includes quasiuniformity. Let , () denote the nodes of lying on the initial smooth surface . The nodes will be evolved in time with the given normal velocity , by solving the ODE
(3) 
which is simply (2) for the nodes . Obviously, the nodes remain on the surface for all times, i.e., for and for all .
Therefore, the smooth surface is also approximated by a discrete surface , whose elements also form a triangulation . We have
As already mentioned in the introduction, we assume that the surface does not develop topological changes due to the evolution. The assumption on quasiuniformity over time, i.e., there is a fixed (independent of ) such that for any triangle the radius of the inscribed circle satisfies
for all , is generally not always satisfied during time evolution.
As an example, from [13], to degenerating surface evolution, we evolved a surface using the ODE (3). As observed in Figure 1: however the initial mesh (left) is quasiuniform, the meshes at later times (middle and right) do not preserve the good mesh qualities. Both quite bad surface resolution and unnecessarily fine elements occur.
3 Computing ALE maps as general constrained problems
We now propose an approach which will be used to determine a suitable ALE map for evolving surfaces, to maintain mesh quality during the surface evolution. In fact, we directly compute the new positions of the nodes.
We consider an other parametrization of , with good mesh properties, and which is different from in (2), called an arbitrary Lagrangian Eulerian map. The corresponding ODE system is
(4) 
with the pure tangential velocity .
For evolving surface problems the surface velocity is usually assumed to be known and normal. However, such pure tangential ALE velocities are not given, and also not easy to obtain, in general. Therefore, we allow velocities (still denoted by) which improve mesh quality, but have small nontangential components, hence may lead points away from the surface. To compensate this, a constraint is introduced in order to keep the surface unaltered. Therefore, the set of ODEs (4) is modified into the following differential algebraic equation (DAE) system (of index 2) with Lagrange multiplier , for ,
(5)  
where . The first argument in indicates that the additional velocity may (and usually does) depend on the whole surface. As it causes no confusion we will drop this argument later on.
The analogous differential algebraic equation system for the nodes of the surface approximation mesh , which are collected into the vector , and with Lagrange multiplier , reads as
(6)  
with initial value . The constraint is meant pointwise, i.e., as for , while the matrix . Concerning notation: we will apply the convention to use boldface letters to denote vectors in or collecting nodal values of discretized variables.
Naturally, the DAE system is independent of the choice of the ALE velocity , if there is some (physical, biological or modeling) knowledge on the general type of the surface evolution a user can propose a suitable ALE velocity accordingly.
In the next sections we will propose a very intuitive way to define the ALE velocity , and we will also discuss numerical methods for the solution of the DAE system.
3.1 A spring system based arbitrary ALE velocity
We use here a simple idea to determine the velocity : let us assume that the nodes of the mesh are connected by springs following the edges of the elements, i.e., the topology of the spring system is determined by the triangulation .
This system defines a force function , which we use to define the ALE velocity by setting
(7) 
The force function is computed based on the connectivity (described by the elements), and by the forces over the edges based on a length function (the desired length of springs). The net force at a node is given by, for ,
(8) 
where the set collects all the edges having as one of their nodes, while simply denotes the other node across the edge , see Figure 2. Then is the force along the edge , given by
and current length .