---

---

1 Introduction↩︎

A wide variety of physical phenomena bring into play evolving interfaces or discontinuities, that take the form of open curves in 2d, or open surfaces in 3d. For instance, in aerodynamics, the wake generated by an aircraft operating at high Reynolds number features vortex sheets, that are infinitely thin open surfaces across which the tangential fluid velocity is discontinuous and whose complex roll-up is driven by self-induced advection [1]. In solid mechanics, fractures are open surfaces across which the displacement of the medium is discontinuous, that evolve under the effect of the stress concentration at their boundaries [2].

A great deal of attention has been paid in the literature to the development of methods for tracking the motion of domains, or equivalently, the closed curves (in 2d) or surfaces (in 3d) defining their boundaries. Broadly speaking, these can be classified into two categories. On the one hand, Lagrangian methods represent the evolving object by a collection of particles and all the attached physical quantities of interest (e.g. the density, velocity and pressure in the case of a fluid, the mass and displacement of a structure) are expressed directly in terms of these, for instance with the help of smoothing kernels. The Smoothed Particle Hydrodynamics (SPH) method, initially proposed in [3] to track bulk domains in the field of astrophysics, then extended to fluid dynamics [4] and solid mechanics [5], is possibly the most famous strategy in this spirit, see [6] for a review. Other Lagrangian methods, such as the Material Point Method [7] or the Immersed Boundary Method [8], complement this representation with a fixed mesh of a larger “hold-all” domain; back-and-forth transfers between the particles and the nodes of the mesh allow to perform mechanical computations on the latter. Lagrangian methods generally suffer from the large number of particles needed to achieve a reasonable accuracy. Furthermore, repeated resampling of the particles is needed to maintain a suitable description of the curve, especially in the regions where the motion induces a large stretching. On the other hand, Eulerian methods capture an evolving object via certain auxiliary fields or markers, that are discretized on a fixed mesh of a computational domain. Among these, the Volume Of Fluid method [9] represents a fluid mixture by the fractions of each fluid inside the grid voxels [10]. Alternatively, the Level Set Method accounts for a moving domain as the negative subdomain of a scalar “level set” function, defined on the whole computational domain, see [11]–[13] and [sec:sec46LSM] below for a brief outline. The Phase-Field Method features a diffuse interface, also captured by a scalar function, whose thickness is penalized by the addition of a perimeter term to the energetic formulation of the evolution [14], [15]. Despite their robust description of large motions, the accuracy of Eulerian representations is limited by the size of the fixed mesh. Furthermore, no explicit support of the interface is available for geometrical or mechanical computations, that have to be expressed in terms of its Eulerian descriptors, often at the expense of approximations. Of course, a whole range of hybrid methods is available, sharing features from the Eulerian and Lagrangian viewpoints. Among these, let us mention the so-called Arbitrary Lagrangian-Eulerian (ALE) methods, where the computational mesh is allowed to move to better track the evolving interface, in a way which may nevertheless differ from the interface motion; see [16] for an overview of this paradigm. Let us also mention the Particle Level Set Method from [17], where a “classical” level set update procedure is helped with a collection of particles evolving in a Lagrangian fashion.

The above investigations primarily address the evolution of closed contours; comparatively little attention has been devoted to the treatment of open objects – that is, curves with endpoints in 2d, as exemplified in 1, or pieces of surfaces having contours in 3d. This task is usually more involved as it demands a careful representation of the motion of the boundary of the object, in addition to that of its bulk structure. Among the Lagrangian strategies implemented to achieve this goal, let us mention the work [18], devoted to tracking the motion of an open curve in two or three space dimensions. The evolving object is discretized with a collection of particles, while a fixed background mesh stores closest-point information at its nodes. The motion of the object is realized by moving the particles, and then updating the closest point information stored at the mesh nodes; a new collection of particles is reconstructed by polynomial fitting, which consistently ensures a uniform sampling of the curve, of the order of the mesh resolution. The information stored at grid points also serves to render topological changes. This framework allows to deal with multiple junctions; however, the treatment of topological changes is heuristic, and no meshed representation of the curve is available. A Eulerian method to track the motion of a 2d open curve is introduced in the work [19], dealing with the numerical simulation of crystal growth: the evolving curve under scrutiny represents the step line of the crystal, and it is driven by a curvature-based velocity field. Elaborating on a suggestion from [20], two level set functions are used to account for the evolving curve, as presented in [sec:sec462LSM]. Such a strategy was used later in [21], [22] for the purpose of image segmentation, and in [23] to fit an open curve (in 2d) or surface (in 3d) to an unorganized point cloud. This idea of using two level set functions for describing an open object has since then been quite popular in the literature, see for instance [24] about its introduction in fracture mechanics, and [25], [26] where it is used to account for the motion of a curve within a surface in 3d. Another Eulerian viewpoint on the evolution of possibly open curves or surfaces is the so-called vector distance function method, introduced in [27]: the considered object is represented by the datum of its unsigned distance function at the vertices of a computational mesh and that of its gradient – which is an extension of the normal vector to its boundary. However elegant, this strategy is a little intricate to implement, as it is difficult to accurately locate the \(0\) level set of the representation. We refer to [28], [29] for implementation details about this elegant viewpoint and [30] for an application to the simulation of crack propagation.

In this open setting also, Eulerian approaches generally do not provide an explicit body-fitted discretization amenable to accurate geometric or mechanical computations, whereas purely Lagrangian approaches struggle with robustness under large deformations and topological changes. The present article aims to reconcile these two requirements by proposing a robust, body-fitted methodology for tracking the motion of open curves in 2d, building on a combination of the Level Set Method and remeshing algorithms. This strategy draws inspiration from our previous work [31]–[33], devoted to the evolution of “bulk” domains, and its extension [34], [35] to the case of regions on surfaces. Our method combines two representations of a two-dimensional open curve \(\Gamma\), or more generally, of a collection of such curves. On the one hand, \(\Gamma\) is meshed explicitly; more precisely, the computational domain \(D\) is equipped with a triangular mesh \({\mathcal{T}}\), and a collection \({\mathcal{L}}\) of edges of \({\mathcal{T}}\) accounts for a discretization of \(\Gamma\). On the other hand, \(\Gamma\) is described implicitly, with the help of two scalar “level set” functions \(\phi, \psi: D \to {\mathbb{R}}\): the \(0\) level set of \(\phi\) is an extension of \(\Gamma\) into a closed curve \(\widetilde{\Gamma}\), and the negative subdomain of \(\psi\) delimits \(\Gamma\) within \(\widetilde{\Gamma}\). Each operation of the evolution workflow (geometric and mechanical computations related to the evaluation of the velocity field, update of the curve...) is applied on the most suitable representation; efficient numerical algorithms allow to pass from one representation to the other whenever needed. This strategy allows to conduct precise mechanical and geometrical calculations related to \(\Gamma\), while leaving room for an arbitrary evolution of the latter.

Most of the conceptual and algorithmic contents of this article hold true in two and three space dimensions, although the implementation is significantly more tedious in the latter case. This article focuses on the 2d setting, where the objects at stake are open curves, or more generally, collections of multiple, disjoint open curves. Although most of the relevant physical applications would take place in 3d, it is already possible to address a few non trivial, interesting applications with this technology. Since the methodology was designed with the three-dimensional setting in mind, the present two-dimensional implementation may be less efficient than more specialized approaches dedicated solely to planar open curves. In a similar spirit, the management of topological changes, however natural with the use of the Level Set Method, is not detailed, as it does not find relevant applications in the presented examples. The extension of the proposed framework to three space dimensions is an ongoing work, that will be presented in the forthcoming article [36]; note however that some of its ingredients have already been used in [37], [38].

The sources code of our numerical implementation are freely available at the following address:

It elaborates on the open-source library developed in [39], dealing with the optimal design of “bulk” shapes.

The remainder of this article is organized as follows. 2 describes the representation of an evolving open curve in \({\mathbb{R}}^2\) by the datum of two level set functions. Then, 3 is devoted to our numerical method for tracking this motion: we outline the key ingredients of our strategy, namely the numerical method for generating two level set functions from a line mesh of an open curve, and, conversely, the remeshing algorithm used to discretize an open curve described by two level set functions. Since these algorithms will be the subject of a dedicated article, in the much more technical three-dimensional setting, we deliberately keep the discussion at a high level, omitting technical details. 4 numerically assesses our framework with two applications where the velocity field at play is “simple”, in the sense that it is explicitly calculated from a discretization of the curve. The next 5 investigate three applications of our methodology in the more challenging perspective of shape optimization, where the velocity field driving the motion typically depends on the calculation of geometrical quantities or on the solution of one or several boundary value problems – thus demanding an accurate, meshed description of the curve. We notably consider the minimization of anisotropic perimeter functionals, the optimization of the laser path of an additive manufacturing process, and the reconstruction of the shape of a fracture within the underground from measurements of the boundary of the upper surface. A conclusion and a few leads for future work are given in 6. This article ends with a series of appendices, containing modeling or technical details related to the different applications of the article; when relevant and not detrimental to clarity, these are provided in a context which is slightly more general than that of the main text.

2 A two-level-set method for the representation of a moving open curve↩︎

This section describes the representation of an evolving open curve by means of two “level set” functions and sets some notations that are used throughout the article. For completeness, we briefly recall in [sec:sec46LSM] the basic features of the “classical” Level Set Method devoted to evolving closed curves, before turning to the two-function variant used in this article to deal with open curves in [sec:sec462LSM].

A brief reminder of the “classical” Level Set Method

The Level Set Method was pioneered in [12] as an efficient means to capture the motion of an evolving domain, and we presently recall its salient features, see for instance the reference books [11], [13] for more exhaustive presentations.

Let \(D \subset {\mathbb{R}}^2\) be a large “hold-all” domain. The Level Set Method consists in capturing a smooth subdomain \(\Omega \subset D\) as the negative region of a scalar “level set” function \(\phi : D \to {\mathbb{R}}\), i.e. \[\label{eq461ls} \forall x \in D, \quad \left\{ \begin{array}{cl} \phi(x) < 0 & \text{if } x \in \Omega,\\ \phi(x) = 0 &\text{if } x \in \partial \Omega, \\ \phi(x) >0 & \text{otherwise.} \end{array} \right.\tag{1}\] This change in perspective does not incur loss of information about \(\Omega\); in particular, all the geometric quantities attached to \(\Omega\) can be expressed in terms of a sufficiently smooth level set function \(\phi\); for instance, the unit normal vector \(n\) to \(\partial \Omega\), pointing outward \(\Omega\), and its mean curvature \(\kappa\) equal: \[n(x) = \frac{\nabla\phi(x)}{\lvert \nabla \phi(x)\lvert} \text{ and } \kappa(x) = \text{\rm div}\left( \frac{\nabla\phi(x)}{\lvert \nabla \phi(x)\lvert}\right), \quad x \in \partial\Omega.\]

Let us now consider a time-dependent domain \(\Omega(t)\), evolving over a time period \((0,T)\) under the effect of a velocity field \(V: {\mathbb{R}}_t \times {\mathbb{R}}^2_x \to {\mathbb{R}}^2\). Introducing a level set function \(\phi(t,\cdot)\) for \(\Omega(t)\), such that 1  holds at each time \(t \in (0,T)\), and a level set function \(\phi_0\) for the initial domain \(\Omega(0)\), an elementary application of the chain rule at the formal level shows that the motion of \(\Omega(t)\) translates as the following advection-like equation for \(\phi\): \[\label{eq461lsadvect} \left\{ \begin{array}{cl} \frac{\partial \phi}{\partial t}(t,x) + V(t,x) \cdot \nabla \phi(t,x) = 0 & \text{on } (0,T) \times D,\\ \phi(0,x) = \phi_0(x) & \text{on } D. \end{array} \right.\tag{2}\] Alternatively, introducing the normal component \(v(t,x) := V(t,x) \cdot \frac{\nabla\phi(t,x)}{\lvert\nabla \phi(t,x)\lvert}\) of the velocity field, 2 takes the form a Hamilton-Jacobi equation: \[\label{eq461lshj} \left\{ \begin{array}{cl} \frac{\partial \phi}{\partial t}(t,x) + v(t,x) \lvert \nabla \phi(t,x) \lvert= 0 & \text{on } (0,T) \times D,\\ \phi(0,x) = \phi_0(x) & \text{on } D. \end{array} \right.\tag{3}\]

Note that 2 3 are not “true” advection or Hamilton-Jacobi equations, insofar as the velocity field \(V(t,x)\) may depend on \(\Omega(t)\), thus on \(\phi(t,x)\) itself. The equations 2 3 hold true in the classical sense as long as \(\Omega(t)\), \(\phi(t,x)\) and \(V(t,x)\) are “smooth”, but they have to be interpreted in the sense of viscosity as soon as singularities appear [40]. We refer to [20], [41] about the mathematical foundations of the level set approach for the mean curvature flow, and to [42] for a larger perspective.

Representation of open curves with two level set functions

This section describes a variant of the Level Set Method adapted to the representation of open curves, whose endpoints are strictly contained in the computational domain \(D\); this idea was proposed in [19], [26], building on a theoretical suggestion of [20], see also [24] about its introduction in the context of fracture mechanics.

Still denoting by \(D\) the fixed computational domain, let \(\Gamma \subset D\) be a smooth, simple open 2d curve. Note that \(\Gamma\) may actually be a collection of several disjoint such curves, but for notational brevity, throughout this article, we refer to it as an “open curve”. We also denote by \(\Sigma\) the boundary of \(\Gamma\), i.e. \(\Sigma\) is a set of isolated endpoints. We characterize \(\Gamma\) via two smooth domains \(\Omega, {\mathcal{O}}\subset {\mathbb{R}}^2\), such that: \[\label{eq462LS} \Gamma = \partial \Omega \cap {\mathcal{O}}, \text{ and } \Sigma = \partial \Omega \cap \partial {\mathcal{O}},\tag{4}\] see 1 for an illustration. We then describe \(\Omega\) and \({\mathcal{O}}\) as the negative subdomains of two respective level set functions \(\phi, \psi: D \to {\mathbb{R}}\), i.e. \[\label{eq462ls} \forall x \in D, \quad \left\{ \begin{array}{cl} \phi(x) < 0 & \text{if } x \in \Omega, \\ \phi(x) = 0 & \text{if } x \in \partial\Omega, \\ \phi(x) > 0 & \text{if } x \in D\setminus \overline{\Omega}, \\ \end{array} \right. \text{ and } \left\{ \begin{array}{cl} \psi(x) < 0 & \text{if } x \in {\mathcal{O}}, \\ \psi(x) = 0 & \text{if } x \in \partial{\mathcal{O}}, \\ \psi(x) > 0 & \text{if } x \in D\setminus \overline{{\mathcal{O}}}. \\ \end{array} \right.\tag{5}\] Intuitively, the \(0\) isoline of \(\phi\) is a closed curve \(\widetilde{\Gamma}:= \partial \Omega\) (or a collection of such), being understood that it may reach the boundary of \(D\) to realize this feature, see 1 (b). The role of the secondary level set function \(\psi\) is to delimit \(\Gamma\) within \(\widetilde{\Gamma}\). In this setting, \(\Gamma\) and \(\Sigma\) are thus characterized by: \[\Gamma = \Big\{ x \in D \text{ s.t. } \phi(x) = 0 \text{ and } \psi(x) <0 \Big\}, \text{ and } \Sigma = \Big\{ x \in D \text{ s.t. } \phi(x) = 0 \text{ and } \psi(x) =0 \Big\}.\]

Analogously to the case of closed curves considered in [sec:sec46LSM], this representation allows to capture the evolution of an open curve \(\Gamma(t)\) over a time period \((0,T)\), according to a velocity \(V: {\mathbb{R}}_t \times {\mathbb{R}}^2_x \to {\mathbb{R}}^2\). Indeed, let \(\phi(t,\cdot)\) and \(\psi(t,\cdot)\) be two level set functions for \(\Gamma(t)\) – i.e. 5 holds at each time \(t\in (0,T)\) – and let \(\phi_0, \psi_0: D \to {\mathbb{R}}\) be two level set functions for the initial set \(\Gamma(0)\). The motion of \(\Gamma(t)\) translates in terms of the following advection-like equations about \(\phi\) and \(\psi\): \[\begin{gather} \label{eq462lsadvect} \left\{ \begin{array}{cl} \frac{\partial \phi}{\partial t}(t,x) + V(t,x) \cdot \nabla \phi(t,x) = 0 & \text{on } (0,T) \times D,\\ \phi(0,x) = \phi_0(x) & \text{on } D, \end{array} \right. \text{ and } \\ \left\{ \begin{array}{cl} \frac{\partial \psi}{\partial t}(t,x) + V(t,x) \cdot \nabla \psi(t,x) = 0 & \text{on } (0,T) \times D,\\ \psi(0,x) = \psi_0(x) & \text{on } D. \end{array} \right. \end{gather}\tag{6}\]

3 Body-fitted evolution of open curves in a two-level set framework↩︎

Let \(D \subset {\mathbb{R}}^2\) be a fixed computational domain, and let \(\Gamma(t) \subset D\) be an open curve, evolving through a time period \((0,T)\) under a velocity field \(V: {\mathbb{R}}_t \times {\mathbb{R}}^2_x \to {\mathbb{R}}^2\). Let \(t^n = n \Delta t\), \(n=0,\ldots, N := T/\Delta t\) be a discretization of \((0,T)\) based on a time step \(\Delta t>0\). For each \(n=0,\ldots,N\), we label with an \(^n\) superscript all the instances of the time-dependent objects under scrutiny at the \(n^{\text{th}}\) iteration.

The proposed strategy for tracking the motion of \(\Gamma(t)\) combines two complementary representations of its configuration \(\Gamma^n\) at each iteration of the evolution process:

• (Meshed representation) The total domain \(D\) is equipped with a valid, conforming and high-quality triangular mesh \({\mathcal{T}}^n\); a discretization of \(\Gamma^n\) is available as a collection \({\mathcal{L}}^n\) of edges pertaining to \({\mathcal{T}}^n\), see 2 (a);

• (Level set representation) \(\Gamma^n\) is accounted for by the datum of two level set functions \(\phi^n, \psi^n :D\to {\mathbb{R}}\), along the lines of [sec:sec462LSM], see 2 (b).

Efficient numerical methods, that are described below (see [sec.gen2LS] [sec.remesh]) allow to pass from one of these representations to the other whenever needed. Thus, each operation of the evolution process can be conducted by using the most suitable representation: the geometrical or mechanical analyses involved in the calculation of the velocity field \(V^n(x)\) enjoy the accuracy of an exact mesh of \(\Gamma^n\), while the evolution of the curve between two successive steps is conveniently accounted for by the Level Set Method.

Our numerical strategy is summarized in 3, and its main steps are discussed in the next sections. In a nutshell, each iteration \(n=0,\ldots\) starts with a meshed representation of \(\Gamma^n\), i.e. \(D\) is equipped with a triangular mesh \({\mathcal{T}}^n\) and \(\Gamma^n\) is discretized as a collection \({\mathcal{L}}^n\) of edges of \({\mathcal{T}}^n\). Two level set functions \(\phi^n\), \(\psi^n\) for \(\Gamma^n\) are then generated at the vertices of the mesh \({\mathcal{T}}^n\), thanks to the algorithm described in [sec:sec46gen2LS]. Then, the velocity field \(V^n(x)\) is calculated; this stage depends on the nature of the motion, and it may require geometric computations (e.g. of the curvature of \(\Gamma^n\)), or the resolution of one or several boundary value problems, which is achieved thanks to the Finite Element Method used on the mesh \({\mathcal{T}}^n\). The update of \(\Gamma^n\) is realized in the level set framework: the equations 6 are solved on \({\mathcal{T}}^n\), with initial data \(\phi^n\), \(\psi^n\), over the time period \((0,\Delta t)\), thanks to the numerical method outlined in [sec:sec46lsadv]. This results in two level set functions \(\phi^{n+1}\) and \(\psi^{n+1}\) for the new curve \(\Gamma^{n+1}\), that are known through their values at the vertices of \({\mathcal{T}}^n\). Eventually, a new mesh \({\mathcal{T}}^{n+1}\) of \(D\) is generated from this datum, where \(\Gamma^{n+1}\) appears explicitly, see [sec:sec46remesh] about this operation.

Generation of two level set functions from a meshed representation

Let the computational domain \(D\) be equipped with a triangular mesh \({\mathcal{T}}\), and let \(\Gamma \subset D\) be an oriented open curve, supplied as a line mesh \({\mathcal{L}}\); for the purpose of this section, the edges of \({\mathcal{L}}\) may not belong to \({\mathcal{T}}\). We aim to generate two level set functions \(\phi\), \(\psi: D \to {\mathbb{R}}\) at the vertices of \({\mathcal{T}}\), representing \(\Gamma\) via 4 and satisfying the orthogonality relation ?? .

This section presents a general and efficient algorithm to achieve this seldom considered task, to the best of our knowledge. Since its more intricate 3d instance will be thoroughly described in a dedicated article [36], we limit ourselves with a brief description in the present 2d context.

The key building block is the celebrated Fast Marching Method [44]–[46], which calculates the (unsigned) distance function \(d(\cdot,K)\) to a subset \(K \subset D\) at the vertices of the mesh \({\mathcal{T}}\) of \(D\). In a nutshell, this method starts with the (unexpensive) calculation of the exact distance to \(K\) at the “close” vertices of the elements of \({\mathcal{T}}\) intersecting \(K\). This “nearby” information is then propagated to the whole mesh \({\mathcal{T}}\): iteratively, the algorithm computes “trial” values from the known distance values and the smallest of them is definitely accepted. Crucially, the method accepts the vertices of \({\mathcal{T}}\) in order, from those closer to those farther from \(K\); this feature has inspired a variant of the method, which concurrently computes the normal extension to the whole mesh \({\mathcal{T}}\) of a quantity defined on \(K\) [47].

Our algorithm for calculating two level set functions \(\phi\) and \(\psi\) representing the open curve \(\Gamma\) proceeds in three stages, that are illustrated on 4.

• Step 1: The unsigned distance function \(d(\cdot,\Gamma)\) is calculated thanks to the Fast Marching Method. In doing so, whenever a vertex \(x \in {\mathcal{T}}\) is accepted by the procedure, the value of the distance is stored in an intermediate level set function \(\phi^{\text{temp}}\), and it is endowed with a sign, depending on the orientation of \(x\) with respect to \(\Gamma\) – an information which is available thanks to the ordered travel of vertices guaranteed by the Fast Marching Method. By construction, the \(0\) level set of \(\phi^{\text{temp}}\) extends \(\Gamma\) into a curve \(\widetilde{\Gamma}\) which is closed, possibly because it goes up to the boundary of the computational domain \(D\), see 1.

• Step 2: The Fast Marching Method is used to calculate the (geodesic) signed distance to \(\Gamma\) within \(\widetilde{\Gamma}\). This information is stored at the vertices of the triangles in the set \({\mathcal{K}}\) of those intersecting the \(0\) level set of \(\phi^{\text{temp}}\). This step yields a function \(\psi^{\text{temp}}\) which contains the values of this signed distance to \(\Gamma\) within \(\widetilde{\Gamma}\), computed at the vertices of the triangles in \({\mathcal{K}}\).

• Step 3: We use once more the Fast Marching Method to compute the signed distance function to \(\widetilde{\Gamma}\), which is stored in \(\phi\). Meanwhile, we extend the function \(\psi^{\text{temp}}\) from \(\widetilde{\Gamma}\) to the whole mesh \({\mathcal{T}}\) in the normal direction to \(\widetilde{\Gamma}\), which yields the secondary level set function \(\psi\).

This strategy guarantees that the curve \(\widetilde{\Gamma}\) extends \(\Gamma\) by a straight line in the neighborhood of \(\Sigma\). Moreover, the orthogonality relation ??  between the \(0\) level sets of \(\phi\) and \(\psi\) is satisfied by construction.

This algorithm is implemented in the open-source library mshdist [48], which is part of the ISCD Toolbox [49].

Resolution of the level set evolution equations

Let \(\Gamma(t)\) be an open curve, evolving through a generic time period \((0,T)\) according to a stationary velocity field \(V(x)\). For any \(t \in (0,T)\), let \(\phi(t,\cdot)\) and \(\psi(t,\cdot)\) be two level set functions for \(\Gamma(t)\), and let \(\phi_0\), \(\psi_0\) be two level set functions for the initial curve \(\Gamma(0)\). We aim to solve the advection equations 6  accounting for the evolutions of \(\phi\) and \(\psi\) on a mesh \({\mathcal{T}}\) of the computational domain \(D\); for definiteness, we focus on that involving \(\phi\): \[\label{eq46advphi} \left\{ \begin{array}{cl} \frac{\partial\phi}{\partial t}(t,x) + V(x) \cdot \nabla \phi(t,x) = 0 & \text{for } t \in (0,T), \: x \in D, \\ \phi(0,x) = \phi_0(x) & \text{for } x \in D, \end{array} \right.\tag{7}\]

In our framework, this equation is solved thanks to the well-known method of characteristics [50], [51]. The latter relies on the following explicit expression of the solution: \[\label{eq46phicharac} \phi(t,x) = \phi_0(X_t(0,x)), \:\: t \in (0,T), \:\: x \in D,\tag{8}\] where, for any time \(t_0\), the mapping \(t\mapsto X_{t_0}(t,x)\) is the characteristic curve of the vector field \(V\) emerging from \(x\) at time \(t_0\), that is: \[\label{eq46characcurveexpl} \left\{ \begin{array}{cl} \frac{\text{\rm d}X_{t_0}}{\text{\rm d}t} (t,x) = V(t,X_{t_0}(t,x)) & \text{for } t \in (0,T), \\ X_{t_0}(t_0,x) = x;& \end{array} \right.\tag{9}\] intuitively, \(X_{t_0}(t,x)\) is the position at time \(t\) of a particle lying in \(x\) at time \(t_0\). Thus, \(\phi\) account for a transport of the quantity \(\phi_0\) along the trajectories induced by the velocity field \(V\).

The numerical resolution of 7 is based on the direct evaluation of the formula 8 at \(t=T\): for each vertex \(x\) of \({\mathcal{T}}\), the ordinary differential equation 9  for the trajectory \(t\mapsto X_T(t,x)\) is solved thanks to a Runge-Kutta 4 method, and the initial function \(\phi_0\) is interpolated at the resulting point from its values at the vertices of \({\mathcal{T}}\).

This algorithm is implemented in the open-source code Advection described in [52], which is part of the ISCD Toolbox [49].

Discretization of an open curve into a mesh from a two-level set representation

Let the computational domain \(D\) be equipped with a triangular mesh \({\mathcal{T}}\), and let \(\Gamma \subset D\) be an open curve. The latter is defined by the datum of two level set functions \(\phi, \psi : D \to {\mathbb{R}}\) that are discretized at the vertices of \({\mathcal{T}}\) and interpolated linearly from these values when evaluated inside the elements of \({\mathcal{T}}\). We aim to create a new, high-quality mesh \(\widetilde{{\mathcal{T}}}\) of \(D\) in which \(\Gamma\) is explicitly discretized.

We proceed in two steps, that are illustrated on 5 and described with a little more details in the next sub-sections.

1. We split each triangle \(T \in {\mathcal{T}}\) intersected by \(\Gamma\) in such a way that the mesh \({\mathcal{T}}^{\text{temp}}\) contains an explicit discretization of \(\Gamma\) as a collection \({\mathcal{L}}^{\text{\rm temp}}\) of edges. This valid and conforming mesh is unfortunately of very low quality.

2. We apply local modifications to \({\mathcal{T}}^{\text{temp}}\) to obtain a new, high-quality mesh \(\widetilde{\mathcal{T}}\) of \(D\) which still features an explicit discretization \(\widetilde{{\mathcal{L}}}\) of \(\Gamma\).

Explicit discretization of an implicitly defined open curve into the computational mesh

This stage is a variation of the marching cubes algorithm [54] and its simplicial version [55], [56], designed for isosurface extraction. It starts from a mesh \({\mathcal{T}}\) of the computational domain \(D\), and the datum of two level set functions \(\phi\), \(\psi\), defined at its vertices, accounting for the curve \(\Gamma\). We form the band \({\mathcal{B}}\) of the elements \(T \in {\mathcal{T}}\) surrounding \(\Gamma\), i.e. satisfying the following two conditions:

(i) A least one vertex of \(T\) has a positive value of \(\phi\), and at least another vertex has a negative value of \(\phi\), i.e. \(T\) intersects the \(0\) level set of \(\phi\);

(ii) At least one vertex of \(T\) bears a negative value of \(\psi\), i.e. \(T\) lies inside or on the boundary of the negative subdomain of \(\psi\).

Then, for each triangle \(T \in {\mathcal{B}}\), we detect the intersection of the isoline \(\left\{ \phi = 0 \right\}\) with the edges of \(T\) by linear interpolation of the values of \(\phi\) at its vertices. A pattern is used to split the triangles sharing at least one of the identified edges by this process in such a way that \(\Gamma\) appears explicitly in the resulting mesh. This operation concerns not only the triangles of \({\mathcal{B}}\), but also the triangles that are adjacent to those in \({\mathcal{B}}\), whose splitting is needed to ensure conformity of the resulting mesh in spite of the fact that they do not intersect the curve \(\Gamma\).

This stage results in a mesh \({\mathcal{T}}^{\text{\rm temp}}\) of \(D\), which is valid, conforming, and contains an explicit discretization of \(\Gamma\) as a collection \({\mathcal{L}}^{\text{\rm temp}}\) of edges. Unfortunately, \({\mathcal{T}}^{\text{\rm temp}}\) is bound to be ill-shaped – i.e. to contain very flat, nearly degenerate elements – since the relative positions of \(\Gamma\) and the vertices of the elements of \({\mathcal{B}}\) are arbitrary, see 5 (b).

Quality-oriented remeshing

In this second stage, we repeatedly apply four local mesh modification operators towards improving the quality of the elements of \({\mathcal{T}}^{\text{\rm temp}}\), while ensuring a fine representation of the geometry of \(\Gamma\):

• Edge split: A “long” edge in the mesh is split into two, after addition of a new vertex in the mesh. All the triangles sharing this edge are split accordingly.

• Edge collapse: One of the two endpoints of a “short” edge is merged with the other, and the vertices of the attached triangles are updated accordingly.

• Edge swap: An edge shared by two triangles is flipped so as to connect the other two vertices of the configuration.

• Vertex relocation: A vertex is moved, while all the connectivities in the mesh remain untouched.

We refer to classical textbooks about meshing such as [57], [58] for a more detailed presentation of these operations.

This stage yields the desired high-quality mesh \(\widetilde{{\mathcal{T}}}\) of \(D\), where a sub-collection \(\widetilde{{\mathcal{L}}}\) of edges explicitly accounts for \(\Gamma\), see 5 (c). Note that this quality-oriented remeshing procedure allows to adapt the size of the elements of \(\widetilde{{\mathcal{T}}}\) to a user-defined size prescription, so as to improve the accuracy of geometric or mechanical computations.

4 Applications examples in physical simulations↩︎

This section presents two “simple” illustrations of our numerical methodology for tracking the motion of an open curve \(\Gamma\) in 2d. All the computations are conducted on a standard Apple MacBookPro laptop with a 2 GHz Quad-Core Intel Core i5 processor and 16 GB of memory. The first [sec:sec46numvor] aims to evaluate its efficiency on the academic problem where \(\Gamma\) evolves under an analytical velocity field. [sec:sec46vortex] then deals with an application in fluid mechanics: \(\Gamma\) represents a two-dimensional vortex sheet whose velocity is given by a weakly singular curve integral.

Evolution of an open curve under an analytical velocity field

This first example aims to appraise the accuracy of our numerical algorithm. It draws inspiration from a classical benchmark test-case proposed in [59] to compare the accuracy of numerical methods for the resolution of the advection equation.

The situation takes places in the 2d unit square \(D = (0,1)^2\). The considered curve \(\Gamma(t)\) evolves through the time period \((0,T)\), with \(T=1.2\), starting from the initial configuration \[\Gamma(0) := \Big\{(0.1+0.8s,0.8) ,\:\: s \in (0,1) \Big\} \cup \Big\{(0.4+0.2s,0.3+0.15s) ,\:\: s\in (0,1)\Big\},\] represented on 6 (a), according to the time-dependent velocity field \(V(t,x)\) defined by: \[\forall t \in (0,T) ,\:\: x= (x_1,x_2) \in D, \quad V(t,x)= \left( \begin{array}{c} \sin(4\pi(x_1+0.5)) \sin(4\pi(x_2-0.3)) \cos(\pi\frac{t}{T}) \\ \cos(4\pi(x_1+0.5)) \cos(4\pi(x_2-0.3)) \cos(\pi\frac{t}{T}) \end{array} \right).\] The latter induces a large deformation of \(\Gamma(t)\) until the time \(t=T/2\); its antisymmetry with respect to \(t=T/2\) then causes the curve to return to its initial configuration at time \(T\).

We apply our numerical 3 to track this motion, choosing a subdivision of \((0,T)\) with \(N= 200\) intermediate times. The remeshing algorithm is required to create elements with minimum and maximum sizes \(0.006\) and \(0.009\), respectively, and the largest mesh produced in the course of the evolution (at \(t=T/2\)) has \(128,729\) vertices. The total computational time equals approximately \(40\) min.

Since the considered motion satisfies \(\Gamma(0)= \Gamma(T)\) at the continuous level, the accuracy of this numerical simulation can be measured in terms of the Hausdorff distance \(d^H(\Gamma(0), \Gamma(T))\) between the initial and final configurations of \(\Gamma(t)\): \[d^H(\Gamma(0), \Gamma(T)) = \max \Big( \rho(\Gamma(0),\Gamma(T)), \rho(\Gamma(T),\Gamma(0)) \Big), \text{ where } \rho(K_1, K_2) := \sup\limits_{x \in K_1} d(x,K_2).\] In the above situation, this error equals \(d^H(\Gamma(0), \Gamma(T)) =1.034\)e\(^{-3}\), which is a much lower value than the minimum size of an edge in the mesh, in spite of the extreme stretching undergone by \(\Gamma\) in the course of the evolution, which validates the accuracy of our method.

Vortex sheet roll-up

In this section, we apply our numerical tracking strategy to the simulation of the dynamics of a vortex sheet in a nearly inviscid 2d fluid; the velocity of such a one-dimensional structure is calculated as a line integral on the evolving curve. The vortex sheet phenomenon has been the subject of extensive investigations: without any claim for exhaustivity, we refer to [60] and Chap. 6 of [61] about its physical modelling, to [60], [62], [63] about its mathematical analysis, and to [64], [65] for numerical simulations by Lagrangian methods; see also [66] for Eulerian simulations of closed vortex sheets via the “classical” Level Set Method.

The phenomenon under scrutiny originally takes place in a 3d medium filled with an incompressible and nearly inviscid fluid. A vortex sheet is a surface along which the velocity “slips”, having discontinuous tangential component and continuous normal component. Such a pattern typically shows up in the wake of an aircraft, or at the limit between two immiscible fluids with different velocities, see 7 (a) for an illustration. In particular configurations, taking advantage of symmetries allows to reduce this situation to that of an open curve \(\Gamma(t)\) evolving within a 2d fluid medium \(D \subset {\mathbb{R}}^2\), as we now consider.

For completeness, a few details about the 2d physical model are sketched in 9. For the purpose of this section, let us solely point out that the dynamics of the fluid and of the discontinuity line \(\Gamma(t)\) are governed by the vortex strength \(\gamma(t,x)\), a scalar quantity defined for \(x\in \Gamma(t)\) which is directly related to the jump in tangential velocity of the fluid. Precisely, the velocity \(u(t,x)\) of the fluid surrounding \(\Gamma(t)\) is given by the following formula: \[\label{eq46velvor} u(t,x) = \frac{1}{2\pi}\int_{\Gamma(t)} \gamma(t,y) \frac{(x-y)^\perp}{\lvert x - y \lvert^2}\:\text{\rm d}s(y), \quad t > 0, \:\: x \in D \setminus \overline{\Gamma(t)}.\tag{10}\] Here, \(v^\perp = (-v_2,v_1)\) stands for the \(90^{°}\) counterclockwise rotate of a vector \(v = (v_1,v_2)\). As throughout the article, we denote by \(\text{\rm d}s\) the integration measure on a codimension \(1\) subset of the plane \({\mathbb{R}}^2\) (i.e. a curve), out of consistency with the general case of a \(d\)-dimensional ambient space.

The velocity field driving the evolution of \(\Gamma(t)\) itself has the same expression 10 , and with a small abuse of notation, it is also denoted by \(u(t,x)\). In this latter case, however, the integrand in 10 is not absolutely integrable, and the formula is understood as a Cauchy principal value. As is customary in the literature, the numerical evaluation of 10 relies on the so-called “vortex-blob” method, which alleviates the singularity of the kernel thanks to the following smooth approximation: \[u(t,x) \approx \frac{1}{2\pi} \int_{\Gamma(t)} \gamma(t,y) \frac{(x-y)^\perp}{\lvert x - y \lvert^2 + \delta^2}\:\text{\rm d}s(y),\] where \(\delta>0\) is a small parameter, see [64].

The description of the evolution of \(\Gamma(t)\) is completed by an equation about the vortex strength \(\gamma(t,x)\). The latter stems from the conservation of circulation within each small portion of \(\Gamma(t)\) as it is driven by the fluid, which reads: \[\lvert \text{\rm com}(\nabla X(t,0,x)) n \lvert \gamma(t,X(t,0,x)) = \gamma(0,x), \quad x \in \Gamma(0),\] where \(t \mapsto X_0(t,x)\) is the characteristic curve of \(u(t,x)\) emerging from a point \(x \in \Gamma(0)\), see 8 .

We apply the methodology of 3 to track the motion of a vortex sheet, in a particular physical configuration considered in [65]. The initial shape of the vortex sheet is the straight line segment \(\Gamma(0) = (-1,1) \times \left\{0\right\}\), and the initial vortex strength is weakly singular at its endpoints: \[\label{eq46inivorstren} \gamma(0,x) = -\frac{x_1}{(1-x_1^2)^{1/2}}, \quad x= (x_1,x_2) \in \Gamma(0).\tag{11}\] The final time of the simulation is \(T=2.2\), and we use the time step \(\Delta t= 0.02\). The vortex-blob parameter is set to \(\delta =0.05\). A few snapshots of the evolution process are represented on 8, which show good agreement with the results of [65]. The maximum number of vertices in a mesh equals \(81,272\), and the total computation takes about \(90\) min.

5 Applications in shape optimization↩︎

This section deals with the numerical resolution of optimization problems where the variable is the shape of an open curve \(\Gamma\). These investigations raise the need for a little background about shape optimization, that we first present in [sec:sec46so]. We then consider three different situations. In [sec:sec46aniper], we optimize the length of an open curve with respect to a Riemannian metric of the plane; the next [sec:sec46am] deals with an application to the 3d printing of a structure: \(\Gamma\) then represents the path of the laser used to raise the temperature of a target region of a metallic powder bed up to the fusion point. Finally, in [sec:sec46fault], we use our methodology to solve an inverse problem of reconstruction of a fault in the underground from observational data.

A primer about shape optimization

The shape optimization problems under scrutiny are of the form \[\label{eq46sopb} \min\limits_{\Gamma \subset D} \: J(\Gamma),\tag{12}\] in which the variable \(\Gamma\) is a 2d open curve (or a collection of such) in a fixed computational domain \(D \subset {\mathbb{R}}^2\). For simplicity, we omit constraints in this formulation, although their presence would not entail any additional conceptual difficulty, see e.g. [67] about this topic.

The treatment of 12 calls for a notion of derivative for a function depending on an open curve, and we rely on the boundary variation method of Hadamard [68]–[72]. In a nutshell, variations of a reference curve \(\Gamma\) are considered of the form: \[\Gamma_\theta = (\text{\rm Id}+ \theta) (\Gamma), \:\: \theta \in {\mathcal{C}}^{1,\infty}({\mathbb{R}}^2;{\mathbb{R}}^2), \:\: \lvert\lvert \theta \lvert\lvert_{{\mathcal{C}}^{1,\infty}({\mathbb{R}}^2;{\mathbb{R}}^2)} < 1,\] i.e. each point \(x\) of \(\Gamma\) is perturbed according to a “small” vector field \(\theta\) in the Banach space \({\mathcal{C}}^{1,\infty}({\mathbb{R}}^2;{\mathbb{R}}^2)\) of bounded vector fields with bounded derivatives. A function \(J(\Gamma)\) is called shape differentiable at \(\Gamma\) if the underlying mapping \(\theta\mapsto J(\Gamma_\theta)\), from \({\mathcal{C}}^{1,\infty}({\mathbb{R}}^2;{\mathbb{R}}^2)\) into \({\mathbb{R}}\), is Fréchet differentiable at \(\theta=0\). This gives rise to the following expansion: \[J(\Gamma_\theta) = J(\Gamma) + J^\prime(\Gamma)(\theta) + \text{\rm o}(\theta) , \text{ where } \frac{\text{\rm o}(\theta)}{\lvert\lvert \theta\lvert\lvert_{{\mathcal{C}}^{1,\infty}(\mathbb{R}^2;\mathbb{R}^2)}} \xrightarrow{\theta\to 0} 0.\] In practice, the knowledge of the shape derivative \(J^\prime(\Gamma)(\theta)\) allows to identify a descent direction for the functional \(J(\Gamma)\), that is, a vector field \(\theta\) such that \(J^\prime(\Gamma)(\theta) < 0\). This property guarantees that for a small enough pseudo-time step \(t>0\) the perturbed configuration \(\Gamma_{t\theta}\) has better performance than \(\Gamma\): \[J(\Gamma_{t\theta}) \:\: \approx \:\: J(\Gamma) + t J^\prime(\Gamma)(\theta) \:\: <\:\: J(\Gamma).\]

The derivatives of the shape functionals considered in this article turn out to be of the form: \[\label{eq46structJp} J^\prime(\Gamma)(\theta) = \int_\Gamma v_\Gamma \theta\cdot n \:\text{\rm d}s + \int_\Sigma \Big( w_\Sigma \theta \cdot n_\Sigma + z_\Sigma \theta \cdot n \Big) \:\text{\rm d}\ell,\tag{13}\] for some known scalar functions \(v_\Gamma : \Gamma \to {\mathbb{R}}\) and \(w_{\Sigma}, z_\Sigma: \Sigma \to {\mathbb{R}}\). This structure reflects that only normal perturbations of the “bulk” of \(\Gamma\) or perturbations of its endpoints may alter the value \(J(\Gamma)\) at first order. A descent direction \(\theta\) can be extracted from 13 by separating the treatments of its tangential and normal components:

• A normal descent direction \(\theta_n := -v_n n\) is obtained thanks to the so-called Hilbertian method, see [68], [73]–[75]: the scalar component \(v_n\) is the solution to the variational problem \[\label{eq46Hilbert} \text{Search for } v_n \in V \text{ s.t. } \forall w \in V, \quad a(v_n,w) = \int_\Gamma v_\Gamma w \:\text{\rm d}s + \int_\Sigma z_\Sigma w \:\text{\rm d}\ell,\tag{14}\] where \(V\) is the Hilbert space \[V = \Big\{ v \in H^1(D) \text{ s.t. } v\lvert_\Gamma \in H^1(\Gamma) \Big\},\] equipped with the following inner product: \[a(v,w) = \alpha^2 \int_D \nabla v \cdot \nabla w \:\text{\rm d}x + \int_D vw \:\text{\rm d}x + \alpha^2 \int_\Gamma \nabla_\Gamma v \cdot \nabla_\Gamma w \:\text{\rm d}s.\] This procedure is justified by the following simple calculation, which directly stems from 13 14 : \[J^\prime(\Gamma)(\theta_n) = -\int_\Gamma v_\Gamma v_n \:\text{\rm d}s - \int_\Sigma z_\Sigma v_n \:\text{\rm d}\ell = -a(v_n,v_n) < 0.\]

• A descent direction \(\theta_\Sigma\) which is tangential to \(\Gamma\) is simply given by: \[\theta_\Sigma = - w_\Sigma n_\Sigma.\]

Combining both observations, the desired descent direction for \(J(\Gamma)\) is then: \[\theta = -v_n n - w_\Sigma n_\Sigma.\]

Note that this separate search for the normal and tangential components of the descent direction \(\theta\) agrees with the representation of \(\Gamma\) by two level set functions \(\phi\) and \(\psi\) described in [sec:sec462LSM]. Indeed, the normal vector field \(\theta_n\) drives the motion of \(\phi\) while leaving \(\psi\) unaltered; on the contrary, advection by the tangential field \(\theta_\Sigma\) does not modify \(\phi\).

Optimization of anisotropic perimeter functionals

This section deals with the minimization of the Riemannian length of a path \(\Gamma\) in the plane. More precisely, we consider the instance of 12 , where the objective function is an anisotropic perimeter functional of the form: \[\label{eq46anisoper} J(\Gamma) = \int_\Gamma \varphi(x,n(x)) \:\text{\rm d}s(x),\tag{15}\] made from a given smooth function \(\varphi :{\mathbb{R}}^2_x \times {\mathbb{R}}^2_n \to {\mathbb{R}}\). The calculation of the shape derivative of such a functional is detailed in 8.

Despite the academic appearance, avatars of this problem show up in multiple applications, such as trajectory optimization [76], car path planning [77] or models for protein folding [78]. Drawing inspiration from [79], we consider an instance of the so-called Zermelo’s navigation problem [80]: a traveler is strolling through a landscape whose topography is defined by a height function \(h:D \to {\mathbb{R}}\) over the 2d unit square \(D\), which is moreover subjected to wind conditions accounted for by the velocity field \(w : D \to {\mathbb{R}}^2\). In this situation, we aim to optimize the trajectory of the traveller, which is represented by an open curve \(\Gamma\) in the domain \(D\), with respect to the length of its image on the graph of \(h\), while taking into account the effort of moving against wind. More precisely, the optimization problem at stake is thus of the form 15 , with \(\varphi\) given by: \[\varphi(x,n) = - n \cdot w^\perp(x) + \sqrt{{\mathcal{M}}(x)n \cdot n},\] where \(w^\perp = (-u_2,u_1)\) is the \(90^{°}\) counterclockwise rotate of \(w\), and \[{\mathcal{M}}(x) = \left(\begin{array}{cc} 1+ \left( \frac{\partial h}{\partial x_2}(x)\right)^2 & - \frac{\partial h}{\partial x_1} (x) \frac{\partial h}{\partial x_2 }(x) \\ -\frac{\partial h}{\partial x_1}(x) \frac{\partial h}{\partial x_2 }(x) & 1+ \left( \frac{\partial h}{\partial x_1} (x)\right)^2 \end{array} \right).\] The first term in the expression of \(\varphi(x,n)\) is minimal when the tangent vector to \(\Gamma\) is aligned with the direction of wind (i.e. \(n\cdot w^\perp\) is maximal) while the second term is the first fundamental form of the landscape expressed in terms of the normal vector to \(\Gamma\).

We consider two different physical situations. In a first example, depicted on 9, the landscape is made of four peaks and no wind is blowing (\(w=0\)). Using a straight line as initial guess, we optimize the path \(\Gamma\) between two fixed endpoints by solving 12 15 . Starting from a straight line, we apply the methodology of 3 to this problem, and the outcome is depicted on 9. The convergence history, reported in 11 (a), shows that the objective function is smoothly minimized until a local minimum of the problem is attained.

In a second experiment, we consider the optimization of the path \(\Gamma\) through a landscape made of three peaks, in the presence of a wind corridor blowing from North-East to South-West: \[w(x) = w_{\text{max}} e^{-\frac{d(x)^2}{\sigma_w^2}} \tau_w , \text{ where } d(x) = \Big(x-(0.1,0)\Big)\cdot n_w \text{ is the distance to the corridor},\] \(w_{\text{max}} = 2\) is the maximum intensity of wind, \(\sigma_w = 0.15\), \(\tau_w = (-0.9923,-0.1240)\) is the direction of wind and \(n_w = \tau_w^\perp\) is the orthogonal direction. The results are depicted on 10, see 11 (b) for the convergence history of the computation. Here, the path takes a detour, occasionally going uphill, in order to avoid moving against the wind.

Optimization of the laser path for additive manufacturing

This section arises in the context of the fabrication of a 3d shape with the Electron Beam Melting 3d printing technology [81], [82]. This process starts by slicing the 3d shape to be produced into a series of horizontal layers, that will be assembled one on top of the other. The construction then takes place in a build chamber which is filled by metallic powder; each layer is assembled by selectively melting the powder in the region occupied by the considered pattern, see 12. Drawing inspiration from the work [83], we aim to optimize the trajectory of the laser during the assembly of a given 2d layer so that powder is melted exactly in the desired region.

We consider the steady-state model proposed in [83], whose physical justification is provided in the appendix of that article. The hold-all domain \(D\) represents the currently processed two-dimensional slice of the build chamber; the laser acts as an instantaneous source of heat along the laser path, which is represented by an open curve \(\Gamma \subset D\). The temperature \(u_\Gamma\) in these circumstances is the \(H^1(D)\) solution to the conductivity equation: \[\label{eq46heateq} \left\{ \begin{array}{cl} -\text{\rm div}(\gamma \nabla u_\Gamma) + \beta (u_\Gamma - u_0) = q \delta_\Gamma& \text{in } D, \\ \gamma \frac{\partial u_\Gamma}{\partial n} = 0 & \text{on } \partial D, \end{array} \right.\tag{16}\] where \(\delta_\Gamma\) stands for integration on \(\Gamma\). Precisely, \(u_\Gamma\) is the solution to the following variational problem: \[\text{Search for } u_\Gamma \in H^1(D)\: \text{ s.t. } \forall v \in H^1(D), \quad \int_D \gamma \nabla u_\Gamma \cdot \nabla v \:\text{\rm d}x + \int_D \beta u_\Gamma v \:\text{\rm d}x = \int_{D} \beta u_0 v \:\text{\rm d}x +\int_\Gamma q v \:\text{\rm d}s.\] In this model, \(\gamma\) is the thermal conductivity of the powder, the parameter \(\beta\) accounts for the heat transmitted from the current 2d layer to the outer medium whose temperature equals \(u_0\), and \(q\) is the power of the laser per unit of length. For simplicity, \(\gamma\), \(\beta\), \(u_0\) and \(q\) are assumed to be constant.

Denoting by \(u_{\text{\rm pc}}\) the phase change temperature of the powder and by \(\Omega_T\) the target shape to be melted within the slice \(D\), we consider the minimization of the following objective function: \[\label{eq46Jam} J(\Gamma) = \int_{\Omega_T}{\left[u_\Gamma - u_{\text{\rm pc}}\right]_-^2 \:\text{\rm d}x} + \int_{D\setminus \overline{\Omega_T}}{\left[u_\Gamma - u_{\text{\rm pc}}\right]_+^2 \:\text{\rm d}x} ,\tag{17}\] where we have set \([t]_+ = \max(0,t)\) and \([t]_- = \min(0,t)\). Intuitively, we wish to optimize the laser path \(\Gamma\) so that the induced temperature should be larger than \(u_{\text{\rm pc}}\) in \(\Omega_T\), thus making the first integral in 17 small, and lower than \(u_{\text{\rm pc}}\) in \(D \setminus \overline{\Omega_T}\), thus making the second integral small.

The shape derivative \(J^\prime(\Gamma)(\theta)\) of this functional is calculated in [83] thanks to the formal Céa’s method [84]; for completeness, the computation is presented in 10 by means of the rigorous technique exposed in e.g. [70], [71].

As such, the minimization of \(J(\Gamma)\) tends to produce curves \(\Gamma\) presenting self-intersections – a pattern that is often undesirable in practice. To alleviate the onset of such features, we follow [85] and bring into play an additional functional: \[C(\Gamma) = \int_\Gamma \Big(\left[d_\Omega(x + d_{\text{\rm min}}n(x))\right]_-^2 + \left[d_\Omega(x - d_{\text{\rm min}}n(x))\right]_+^2 \Big)\:\text{\rm d}s,\] where \(\Omega\) denotes any smooth bounded domain whose boundary encloses \(\Gamma\) and \(d_\Omega\) stands for the signed distance function to \(\Omega\), see ?? . Loosely speaking, the minimization of \(C(\Gamma)\) imposes that for each parameter value \(t \in (-\frac{d_{\text{\rm min}}}{2},\frac{d_{\text{\rm min}}}{2})\), the offset curve \(\Gamma_{tn} = (\text{\rm Id}+ tn)(\Gamma)\) should be diffeomorphic to \(\Gamma\): \(C(\Gamma)\) takes large values if at some point \(x \in \Gamma\), one of the normal rays \((x,x + d_{\text{\rm min}}n(x))\) or \((x,x - d_{\text{\rm min}}n(x))\) crosses the boundary of \(\Omega\). The calculation of the shape derivative of \(C(\Gamma)\) and its numerical integration into a shape optimization algorithm are achieved exactly along the lines of our previous work [86], which is devoted to a variational method for enforcing geometric constraints, expressed in terms of the signed distance function to the optimized shape.

All things considered, the shape optimization problem under scrutiny reads: \[\label{eq46sopbam} \min\limits_{\Gamma} \:\: \Big( J(\Gamma) + \ell C(\Gamma) \Big),\tag{18}\] where \(\ell\) is a fixed weight.

We apply this shape optimization framework to the setting exemplified in 13 (a): the design domain \(D\) is the unit square \((0,1)^2\) and the shape \(\Omega_T\) of the 2d slice to be assembled is a slightly smaller square. Starting from the simple initial curve \(\Gamma^0\) depicted on 14 (a), we apply our numerical 3 to the resolution of the optimization problem 18 . The parameters of the computation are \(\gamma =0.01\), \(\beta = 0.5\), \(q = 100\), \(u_0 = 0\), \(u_{\text{\rm pc}}= 350\) and \(\ell =1\text{e}9\). The maximum number of vertices in a mesh is \(48,245\) (for about twice as many elements) and the total computational time is about 20 mn. A few snapshots of the optimization process are reported on 14 and the convergence history in 13 (b) shows the smooth decrease of the objective function until a local minimum is attained.

We next turn to a second example, associated to the more complex layer pattern depicted on 15 (a), using the same physical parameters as for the previous example; a few intermediate configurations of the optimized curve are presented on 16 and the convergence history of the computation is presented on 15 (b).  

Least-square reconstruction of a fault line in the underground

This section deals with an application of our framework in geophysics. The considered open curve \(\Gamma\) accounts for a fault line in the underground Earth, where the displacement is discontinuous. We aim to reconstruct the shape of \(\Gamma\) by minimizing a least-square function of the difference between the displacement predicted by a “forward” mathematical model and that observed in the course of a geological event.

We first present the exact mathematical model predicting the displacement of the underground from the knowledge of the fracture \(\Gamma\) in [sec:sec46modfault], before introducing an approximate version in [sec:sec46faultdirect] which simplifies the numerical resolution. The inverse reconstruction problem of \(\Gamma\) from observations of the displacement is eventually addressed in [sec:sec46IPfracture].

Description of the physical model

Let \(D\subset {\mathbb{R}}^2\) be a “hold-all” domain, accounting for a two-dimensional section of the geographical region of interest, which is assumed to be made of a linearly elastic material. This region contains a fault, modeled as an open curve \(\Gamma \Subset D\), whose unknown shape is to be determined, see 17 (a). Along the fault, the elastic displacement is discontinuous (i.e. the underground may slip differently on both sides), but the normal traction is continuous. The bottom side \(\partial D_D\) of \(D\) is assumed to be fixed and the remaining boundary \(\partial D \setminus \overline{\partial D_D}\), which is composed of the upper surface \({\Gamma_{\text{\rm obs}}}\) and the lateral sides, is traction-free, i.e. homogeneous Neumann boundary conditions are applied. The displacement of the crust is the unique solution \(u_\Gamma \in H^1(D \setminus \Gamma)^2\) to the linear elasticity system: \[\label{eq46elasfrac} \left\{ \begin{array}{cl} -\text{\rm div}(Ae(u_\Gamma)) = 0 & \text{in } D \setminus \overline{\Gamma}, \\ u_\Gamma = 0 &\text{on } \partial D_D, \\ Ae(u_\Gamma)n = 0 & \text{on } \partial D \setminus \overline{\partial D_D},\\ \left[u_\Gamma\right] = g_\Gamma \text{ and } \left[Ae(u_\Gamma)n\right] = 0 & \text{on } \Gamma. \end{array} \right.\tag{19}\] Here, \(e(u) = \frac{1}{2}(\nabla u+ \nabla u^T)\) is the strain tensor associated to a displacement field \(u : D \to {\mathbb{R}}^2\) and \(A\) is the Hooke’s law of the elastic material in the crust: \[\text{For all } 2\times 2 \text{ symmetric matrix } e, \quad Ae = 2\mu e + \lambda \text{\rm tr}(e) \text{\rm I},\] where \(\lambda = 0.5769\) and \(\mu = 0.3846\) are the Lamé coefficients. In 19 , \(n\) stands for the unit normal vector to \(\Gamma\); we also denote by: \[\label{eq46defjump} \alpha^+(x) = \lim\limits_{t \to 0 \atop t > 0} \alpha(x+tn(x)), \:\: \alpha^-(x) = \lim\limits_{t \to 0 \atop t > 0} \alpha(x-tn(x)), \text{ and } \left[ \alpha \right] (x) = \alpha^+(x) - \alpha^-(x)\tag{20}\] the one-sided limits and the jump of a quantity \(\alpha\) which is smooth enough from either side of \(\Gamma\). The so-called slip vector \(g_\Gamma\) belongs to the functional space \[\label{eq46H12t} \widetilde{H}^{1/2}(\Gamma)^2 := \left\{ u \in L^2(\Gamma)^2 \text{ s.t. } \widetilde{u} \in H^1(\widetilde{\Gamma})^2 \right\},\tag{21}\] where \(\widetilde{\Gamma}\) is an arbitrary extension of \(\Gamma\) into a closed curve and \(\widetilde{u}\) is the extension of \(u\) to \(\widetilde{\Gamma}\) by \(0\). The expression of \(g_\Gamma\) depends on the shape of \(\Gamma\). Let \(c_0, c_1\) be the endpoints of \(\Gamma\) and, for any points \(x_0, x_1 \in \Gamma\), let \(\Gamma_{x_0,x_1}\) denote the region of \(\Gamma\) comprised between \(x_0\) and \(x_1\); \(g_\Gamma(x)\) is of the form: \[\label{eq46slipvec} \forall x \in \Gamma, \quad g_\Gamma(x) = g(s_\Gamma(x)),\tag{22}\] where \(g : [0,1] \to {\mathbb{R}}^2\) is a given smooth function vanishing at \(0\) and \(1\), and \[s_\Gamma(x) := \frac{\lvert \Gamma_{c_0,x}\lvert}{\lvert\Gamma\lvert}, \text{ where } \lvert \Gamma \lvert = \int_\Gamma \:\text{\rm d}s,\] is the normalized arc length of the region \(\Gamma_{c_0,x}\) of \(\Gamma\) delimited by \(c_0\) and \(x\).

Numerical simulation

The numerical simulation of 19 is cumbersome, even when a body-fitted description of the fracture \(\Gamma\) is available as a collection \({\mathcal{L}}\) of edges of the mesh \({\mathcal{T}}\) of the computational domain \(D\). Indeed, “classical” continuous finite element methods cannot easily represent discontinuous functions: one would have to duplicate the nodes on \(\Gamma\) into different copies, associated to the triangles located on either side of the latter, and then to modifiy the element-to-vertex relations accordingly. This procedure is quite intrusive with respect to the implementation of the finite element solver, see nevertheless [38] for a numerical algorithm allowing to do so.

Here, we rather rely on the elegant method introduced in [88], [89], based on a duplication of the sought function. Briefly, let us introduce two subdomains \(\Omega_1\), \(\Omega_2 \subset D\) dividing \(D\) according to \(\Gamma\), that is: \[\Omega_1 \cap \Omega_2 = \emptyset , \:\: \overline{D} = \overline{\Omega_1} \cup \overline{\Omega_2}, \text{ and } \Gamma \subset \partial \Omega_1\cap \partial \Omega_2.\] The subdomains \(\Omega_1\) and \(\Omega_2\) lie respectively “below” and “above” the fracture set \(\Gamma\): the normal vector \(n\) is oriented from \(\Omega_1\) to \(\Omega_2\), see 17 (a). For instance, in the numerical context where \(\Gamma\) is accounted for by two level set functions \(\phi,\psi : D \to {\mathbb{R}}\) as in 4 , \(\Omega_1\) and \(\Omega_2\) may be defined as the negative and positive subdomains of \(\phi\), respectively.

Let us then introduce the functional space \[H^1_{\partial D_D}(D) = \left\{ u \in H^1(D) \text{ s.t. } u = 0 \text{ on } \partial D_D \right\},\] and consider the following variational problem: \[\begin{gather} \label{eq46appvarffract} \text{Search for } (u_{1,\varepsilon}, u_{2,\varepsilon}) \in H^1_{\partial D_D}(D)^4 \text{ s.t. } \forall (v_1, v_2) \in H^1_{\partial D_D}(D)^4, \\ \int_D A_{1,\varepsilon} e(u_{1,\varepsilon}) : e( v_1) \:\text{\rm d}x + \int_D A_{2,\varepsilon} e(u_{2,\varepsilon}) : e(v_2) \:\text{\rm d}x + \frac{1}{\varepsilon} \int_\Gamma ( u_{2,\varepsilon} - u_{1,\varepsilon} ) \cdot (v_2-v_1) \:\text{\rm d}s =\\ \frac{1}{\varepsilon} \int_\Gamma g_\Gamma \cdot (v_2-v_1) \:\text{\rm d}s, \end{gather}\tag{23}\] where we have defined the extensions of the Hooke’s law \(A\) outside \(\Omega_1\) and \(\Omega_2\) by: \[A_{i,\varepsilon}(x) = \left\{ \begin{array}{cl} A & \text{if } x \in \Omega_i, \\ \varepsilon A & \text{otherwise}. \end{array} \right.\] Loosely speaking, 23 is obtained from 19 by duplicating the sought displacement \(u\) (and the test function \(v\)) into the pair \((u_{1,\varepsilon}, u_{2,\varepsilon})\), where \(u_{i,\varepsilon}\) is meant as an approximation of the restriction of \(u\) to \(\Omega_i\), \(i=1,2\). Both \(u_{1,\varepsilon}\) and \(u_{2,\varepsilon}\) are continuous on \(D\) but only their respective values on \(\Omega_1\) and \(\Omega_2\) are relevant. Furthermore, their traces do not agree on \(\Gamma\): the difference \(u_{2,\varepsilon} - u_{1,\varepsilon}\) is imposed to be an approximation of \(g_\Gamma\) by penalization. The justification of this approximation is relatively classical; for completeness, the argument is sketched in [sec:app46justifbroken], in the model setting of the conductivity equation, as a variation of the analysis conducted in [88].

An application example of this approximation strategy is presented on 17 (b), where the solution \((u_{1,\varepsilon}, u_{2,\varepsilon})\) to 23 is computed in the physical situation of 17 (a). The penalization parameter is set to \(\varepsilon= 1\text{e}^{-15}\) and the slip profile \(g\) is defined by: \[\forall s \in [0,1], \quad g(s) = \left\{ \begin{array}{cl} \left(0.02\left(1-\left(\frac{s-0.5}{0.49}\right)^4 \right),0.05 \left(1-\left(\frac{s-0.5}{0.49}\right)^4 \right) \right) & \text{if } s \in [0.01,0.99],\\ 0 & \text{otherwise} \end{array} \right.\]

Reconstruction of the fracture \(\Gamma\) from observational data

This section deals with the inverse problem associated to the “forward” model of [sec.modfault] [sec.faultdirect] for the elastic displacement of the underground in the presence of a fracture. We aim to reconstruct the shape of a fracture \(\Gamma\) within the 2d elastic medium \(D\) from the observation of the displacement \(u_{\text{\rm obs}}: {\Gamma_{\text{\rm obs}}}\to {\mathbb{R}}^2\) of the surface \({\Gamma_{\text{\rm obs}}}\).

We follow the approach in [90]. The fracture \(\Gamma\) is sought as a minimizer of the average on \({\Gamma_{\text{\rm obs}}}\) of the least-square difference between the observed displacement \(u_{\text{\rm obs}}\) and the model prediction \(u_\Gamma\), solution to 19 . Precisely, we consider the optimization problem: \[\label{eq46sopbfrac} \min\limits_{\Gamma \subset D} J(\Gamma), \text{ where } J(\Gamma) = \frac{1}{2} \int_{{\Gamma_{\text{\rm obs}}}} \lvert u_\Gamma -u_{\text{\rm obs}}\lvert^2 \:\text{\rm d}s.\tag{24}\] The shape derivative of \(J(\Gamma)\) is calculated in 7.

The physical situation at stake is depicted on 18 (a); the observed displacement \(u_{\text{\rm obs}}: {\Gamma_{\text{\rm obs}}}\to {\mathbb{R}}^2\) of the surface is that caused by the known fracture pattern \(\Gamma_T\) in there. For a given shape \(\Gamma\) of the fracture, the displacement \(u_\Gamma\) is approximated via the couple \((u_{1,\varepsilon},u_{2,\varepsilon})\), solution to 23 , along the lines of [sec:sec46faultdirect]. The reconstruction problem 24 is known to be very instable, especially when the slip vector vanishes at the endpoints of the fracture and only one measurement is used, see the discussion in [90] where this phenomenon is observed even though the reconstructed fractured is constrained to be a line segment. Here, we content ourselves with observing that the reconstructed fracture roughly has a similar inclination as the exact one, and that the objective function \(J(\Gamma)\) is dutifully minimized: its value is decreased by 500 times within 100 iterations.

6 Conclusion and perspectives↩︎

In this article, we have proposed a novel numerical framework for tracking the motion of an evolving open curve, or a collection of such. This strategy combines two representations of the curve at each iteration of the evolution process: (i) an explicit discretization, as a subset of the edges of the computational mesh, which allows to perform accurate evaluations of its geometric features and precise finite element computations associated to its physical behavior; (ii) an implicit representation, via a variant of the Level Set Method, which allows to capture its evolution in a robust and efficient manner, however dramatic. We have used our framework to address various problems, such as the simulation of a vortex sheet roll-up, the optimization of the path of the laser ensuring the 3d printing of a shape by the Electron Beam Melting technology, or the reconstruction of a fault in the underground from observational data.

This work is preliminary in many respects and calls for multiple perspectives. At first, the proposed methodology can be applied to various physical problems, beyond those exemplified in this article:

• Open curves are ubiquitous in image segmentation, where they can be optimized to best indicate sharp variations of intensity or color, see e.g. [18], [91].

• A similar formalism to that of the path optimization problem discussed in [sec:sec46am] could be applied to the optimal graft of a thin ligament to a structure in order to improve its robustness: our previous works [92], [93] indeed reduce this task to the minimization of an anisotropic perimeter functional after calculation of a “topological ligament expansion” of a mechanical shape functional of interest.

• We aim to apply the methodology in the physical setting of electromagnetism, to optimize the placement of current lines within the ferromagnetic region of an electric motor.

Another lead for improving the proposed framework concerns the handling of one-dimensional structures showing branching or multiple junctions: its coupling with a strategy for decomposing a branched structure into a collection of manifold curves would allow to simulate complex fracture propagation patterns.

Finally, and perhaps most importantly, as we have mentioned in the introduction, the original motivation and natural perspective of this work is to develop an efficient and accurate, body-fitted tracking strategy for the motion of open surfaces in three space dimensions, with applications including the optimization of shells or electromagnetic screens. Although the implementation effort is then expected to be more intense, the methodology of this article was developed with this ambition in mind, and its extension to this setting does not present any theoretical obstruction.

Acknowledgements. This work is partially supported by the projects ANR-24-CE40-2216 STOIQUES, ANR-22-CE46-0006 StableProxies and CNRS-MITI DOLMEN. Part of it was realized while the author was visiting the Department of Mathematics of the University of Trento, whose hospitality is gratefully acknowledged.

7 A few facts from tangential calculus↩︎

This appendix gathers some basic facts from differential calculus on a smooth hypersurface of \({\mathbb{R}}^d\), \(d \geq 2\). We mainly follow [70], see also classical books such as [94], [95] for differential calculus on manifolds.

Let \(\Gamma\) be an oriented open surface in \({\mathbb{R}}^d\), \(d \geq 2\), whose contour \(\Sigma := \partial \Gamma\) is then a \((d-2)\)-dimensional closed submanifold of \({\mathbb{R}}^d\). We denote by \(n : \Gamma \to {\mathbb{R}}^d\) the unit normal vector, whose orientation is consistent with that of \(\Gamma\). Moreover,

• The tangential gradient of a function \(u : \Gamma \to {\mathbb{R}}\) of class \({\mathcal{C}}^1\) is defined by \(\nabla_\Gamma u := \nabla \widetilde{u} - (\nabla \widetilde{u} \cdot n ) n\), where \(\widetilde{u}\) is an arbitrary extension of \(u\) to an open neighborhood of \(\Gamma\) in \({\mathbb{R}}^d\).

• The tangential divergence of a vector field \(v : \Gamma \to {\mathbb{R}}^d\) is defined by \(\text{\rm div}_\Gamma(v) = \text{\rm div}(\widetilde{v}) - \nabla \widetilde{v} n \cdot n\), where \(\widetilde{v}\) is an arbitrary extension of \(v\) to an open neighborhood of \(\Gamma\) in \({\mathbb{R}}^d\).

Let us now recall the following formula for changing variables in surface integrals.

We now state a well-known integration by parts formula on a surface.

8 Calculation of the shape derivative of anisotropic perimeter functionals↩︎

This appendix details the calculation of the shape derivatives of the anisotropic perimeter functionals considered in [sec:sec46aniper]. Placing ourselves in the general context of a \(d\)-dimensional ambient space, we consider a functional of the form \[J(\Gamma ) =\int_\Gamma \varphi(x, n_\Gamma(x)) \:\text{\rm d}s,\] depending on an open smooth hypersurface \(\Gamma\) through its unit normal vector \(n_\Gamma : \Gamma \to {\mathbb{R}}^d\), that we simply denote by \(n\) when the concerned hypersurface is clear. Here, \(\varphi : {\mathbb{R}}^d_x \times {\mathbb{R}}^d_n \times {\mathbb{R}}\) is a given smooth function; the gradients of the partial mappings \(x \mapsto \varphi(x,n)\) and \(n \mapsto \varphi(x,n)\) are respectively denoted by \(\nabla_x \varphi (x,n)\) and \(\nabla_n \varphi(x,n)\).

 

9 Additional details about the 2d vortex sheet evolution problem↩︎

This appendix expands a little on the vortex sheet roll-up example considered in [sec:sec46vortex] and its numerical treatment. The mathematical model is described in a formal way in [sec:app46vormod] and a few additional details about its practical implementation are provided in [sec:app46vorimp].

Mathematical model of the vortex sheet evolution in 2d

A vortex sheet is a phenomenon that may develop within an incompressible and nearly inviscid fluid, i.e. when the Reynolds number is very high. As we have mentioned in [sec:sec46vortex], it takes the form of a curve (in 2d) or a surface (in 3d) \(\Gamma(t)\) across which the fluid velocity “slips”, i.e. its normal component is continuous, but its tangential component is discontinuous. Such a structure has its own dynamics: it expands and rolls up under the effect of the complex velocity patterns that it itself creates. We refer to [1], [96] for general references about fluid mechanics, to [97], [98] about its mathematical framework, and to [60]–[63] about the specific subject of vortex patterns.

Formally, let us assume that the fluid under scrutiny is contained in the infinite plane \({\mathbb{R}}^2\), i.e. we neglect the boundary effects induced by the use of a fixed computational domain \(D\). We denote by \(u : {\mathbb{R}}_t \times {\mathbb{R}}_x^2 \to {\mathbb{R}}^2\) and \(p:{\mathbb{R}}_t \times {\mathbb{R}}_x^2 \to {\mathbb{R}}\) the fluid velocity and pressure, respectively. Denoting by \(T\) the final time of the study, the couple \((u,p)\) is solution to the inviscid, incompressible Navier-Stokes equations: \[\label{eq46NS} \left\{ \begin{array}{cl} \frac{\partial u}{\partial t} + \nabla u u + \nabla p = 0 & \text{for } (t,x) \in (0,T) \times {\mathbb{R}}^2, \\ \text{\rm div}(u) = 0 & \text{for } (t,x) \in (0,T) \times {\mathbb{R}}^2, \\ u(0,x) = u_0(x) & \text{for } x \in {\mathbb{R}}^2. \end{array} \right.\tag{29}\] The first equation in 29 accounts for the balance of momentum within the fluid, while the second one accounts for its incompressibility; the initial velocity profile \(u_0: {\mathbb{R}}^2 \to {\mathbb{R}}^2\) is assumed to be known.

A handful quantity in the description of this situation is the vorticity of the fluid \(\omega : {\mathbb{R}}_t \times {\mathbb{R}}_x^2 \to {\mathbb{R}}\), which is defined by: \[\label{eq46omcurl} \omega(t,x) = \text{\rm curl}(u)(t,x) := \frac{\partial u_2}{\partial x_1}(t,x) - \frac{\partial u_1}{\partial x_2} (t,x).\tag{30}\] Intuitively, at each time \(t\), \(\omega(t,x)\) accounts for the infinitesimal rotating motion of the fluid around \(x\). By taking the curl of the balance of momentum relation and using the incompressibility condition \(\text{\rm div}(u) = 0\), an elementary calculation yields the following equation about \(\omega\): \[\frac{\partial \omega}{\partial t} + u\cdot \nabla \omega =0,\] which expresses that the vorticity is transported along with the fluid particles.

Another important quantity in the description of vortex sheets is the stream function \(\psi: {\mathbb{R}}_t \times {\mathbb{R}}_x^2 \to {\mathbb{R}}\) of the fluid. Since \(\text{\rm div}(u) = 0\), the Helmholtz decomposition implies that there exists a function \(\psi\), which is unique up to a constant, such that: \[\label{eq46defpsi} u = \nabla^\perp \psi := \left(-\frac{\partial \psi}{\partial x_2},\frac{\partial\psi}{\partial x_1}\right).\tag{31}\] Combining this identity with the definition 30 of \(\omega\), we arrive at: \[\label{eq46Deltapsi} \Delta \psi = \omega.\tag{32}\]

Let us now assume that a vortex sheet is present within the fluid, that takes the form of an oriented open curve \(\Gamma(t)\). We denote by \(n_t\) the unit normal vector to \(\Gamma(t)\) and by \(\tau_t\) its unit tangent vector, such that \((\tau_t,n_t)\) is a direct orthonormal frame of the plane, i.e. \(n_t = \tau_t^\perp := (-\tau_{t,2,}\tau_{t,1})\), see 7 (b). Let us recall from [sec:sec46fault] that, when \(\alpha\) is a quantity which is smooth on either side of \(\Gamma(t)\), but possibly discontinuous across \(\Gamma(t)\), we denote by \[\alpha^\pm(x) = \lim\limits_{s \to 0 \atop s > 0} \alpha(x \pm s n_t(x)) \text{ and } \left[ \alpha \right] (x) = \alpha^+(x) - \alpha^-(x)\] the one-sided limits and the jump of \(\alpha\) at \(x \in \Gamma(t)\), respectively. Mathematically, the vortex sheet pattern is characterized by the fact that, at each time \(t\), the vorticity is concentrated on \(\Gamma(t)\), i.e. \[\omega(t,x) = \gamma(t,x) \delta_{\Gamma(t)},\] where \(\delta_{\Gamma}\) is the distribution accounting for integration over \(\Gamma(t)\), and the scalar factor \(\gamma(t,x)\) is called the strength of the vortex. This relation has to be understood in the sense of distributions in \({\mathbb{R}}^2\): for any open subset \(U \subset {\mathbb{R}}^2\), it holds: \[\label{eq46locvor} \int_U \omega(t,x) \:\text{\rm d}x = \int_{\Gamma(t) \cap U} \gamma(t,x) \:\text{\rm d}s.\tag{33}\] Substituting 32 for \(\omega(t,x)\) in 33  and integrating by parts, we arrive at the following boundary-value problem for the stream function \(\psi(t,\cdot)\) at each time \(t>0\): \[\label{eq46bvppsi} \left\{ \begin{array}{cl} \Delta \psi = 0 & \text{in } {\mathbb{R}}^2 \setminus \Gamma(t), \\ -\left[ \frac{\partial \psi}{\partial n}\right] = \gamma(t,\cdot) & \text{on } \Gamma(t), \end{array} \right.\tag{34}\] where the sign comes from the chosen orientation in the definition of the jump across \(\Gamma(t)\). It follows from 34  that \(\psi(t,\cdot)\) can be expressed as a single layer potential, \[\psi(t,x) = -\int_{\Gamma(t)} \gamma(t,y) G(x,y) \:\text{\rm d}s(y),\] where \(G(x,y) = \frac{1}{2\pi} \log\lvert x - y \lvert\) is the fundamental solution of the Laplace operator in 2d, see for instance [99], [100] about the topic of potential theory. Using once again the relation 31 between \(u\) and \(\psi\), we arrive at the following representation formula for \(u\) in terms of \(\gamma\): \[u(t,x) = \int_{\Gamma(t)} \gamma(t,y) K(x,y) \:\text{\rm d}s(y), \text{ where } K(x,y) := -\frac{1}{2\pi} \frac{ (x-y)^\perp}{\lvert x - y\lvert^2}.\] The classical jump relations for the single layer potential readily imply that \(u\) has continuous normal component across \(\Gamma(t)\), but that its tangential component is discontinuous: \[\label{eq46jump} \left[u\cdot n_t \right] = 0 ,\text{ and } (u \cdot \tau_t)^\pm= \pm \frac{1}{2} \gamma(t,x) + \text{\rm p.v.}\int_{\Gamma(t)} \gamma(t,y) K(x,y) \cdot \tau_t(x) \:\text{\rm d}s(y),\tag{35}\] where the last integral is understood in the sense of a Cauchy principal value, see again [99], [100]. In particular, this shows that the vortex strength \(\gamma(t,x)\) coincides with the jump \(-[u\cdot \tau_t]\) in tangential velocity.

Let us now describe the evolution of \(\Gamma(t)\). It can be shown that the consistency of the latter with the inviscid Navier-Stokes equations 29 for the surrounding fluid imposes that the velocity field driving its motion should be equal to the average of the one-sided values of the fluid velocity, see Chap. 6 of [61]. Using the jump relations 35 , this means that the velocity of \(\Gamma(t)\), that we still denote by \(u(t,x)\) with a small abuse of notation, is given by the so-called Birkhoff-Rott equation: \[\label{eq46pvuvortex} u(t,x) = \text{\rm p.v.}\int_{\Gamma(t)} \gamma(t,y) K(x,y) \:\text{\rm d}s(y).\tag{36}\] To complete this model, we need to characterize the evolution of the vortex strength \(\gamma(t,x)\). This task involves another conservation principle: the total vortex strength contained in any portion of \(\Gamma(t)\) is conserved during the motion. Mathematically, for any \(s >0\), let us denote by \(t\mapsto X_s(t,x)\) the characteristic curve of the velocity field \(u(t,x)\) emerging from \(x\) at time \(s\), which is the unique solution to the following ordinary differential equation, see 9 : \[\label{eq46odecharacbis} \left\{ \begin{array}{cl} \frac{\text{\rm d}X_s}{\text{\rm d}t}(t,x) = u(t,X_s(t,x)) & \text{for } t \in (0,T), \\ X_s(s,x) = x.& \end{array} \right.\tag{37}\] Mathematically, the above conservation principle means that, for any region \(G \subset \Gamma(0)\), it holds: \[\frac{\text{\rm d}}{\text{\rm d}t} \left( \int_{X_0(t,G)} \gamma(t,y) \:\text{\rm d}s(y) \right) = 0.\] Using the change of variables in surface integrals recalled in 1, this rewrites: \[\frac{\text{\rm d}}{\text{\rm d}t} \left( \int_{G} \lvert \text{\rm com}\nabla X_0(t,y) n_t(y) \lvert \gamma(t,X_0(t,y)) \:\text{\rm d}s(y) \right) = 0.\] Finally, since this property holds for any subset \(G \subset \Gamma(0)\), the quantity \(\lvert \text{\rm com}\nabla X_0(t,y) n_t(y) \lvert \gamma(t,X_0(t,y))\) is conserved, and so: \[\label{eq46circconserv} \text{For all } x \in \Gamma(0), \quad \lvert \nabla X_0(t,x) \tau_0(x) \lvert \gamma(t,X_0(t,x)) = \gamma(0,x),\tag{38}\] where we have used the following elementary calculation: \[\lvert \text{\rm com}\nabla X_0(t,x) n_0(x) \lvert = \lvert \nabla X_0(t,x) \tau_0(x) \lvert.\] Equivalently, 38  has the following equivalent backward expression: \[\label{eq46circconservback} \text{For all } y \in \Gamma(t), \quad \gamma(t,y) = \frac{\gamma(0,X_t(0,y))}{\lvert \nabla X_0(t,X_t(0,y)) \tau_0(X_t(0,y))\lvert} = \frac{\gamma(0,X_t(0,y))}{\lvert \nabla X_t(0,y)^{-1} \tau_0(X_t(0,y))\lvert}.\tag{39}\] Finally, the equations 36 38 (or 39 ) completely characterize the dynamics of the vortex sheet \(\Gamma(t)\).

A few implementation details

The numerical implementation of the above model for the simulation of vortex sheets raises a few practical issues, that we broach in this section.

Update of the vortex strength

At each iteration \(n=0,\ldots\) of the evolution, the vortex strength \(\gamma(t^n,x)\) is discretized at the vertices of the line mesh \({\mathcal{L}}^n\) of \(\Gamma(t^n)\), whose edges explicitly appear in the mesh \({\mathcal{T}}^n\) of \(D\), as one key feature of our body-fitted tracking method, see 3. The calculation of the updated values \(\gamma(t^{n+1},\cdot)\) at the vertices of the mesh \({\mathcal{L}}^{n+1}\) for \(\Gamma(t^{n+1})\) from those of \(\gamma(t^n,\cdot)\) on \(\Gamma(t^n)\) is based on the formula 39 , which can be given an iterative flavor: \[\label{eq46vorstrenupdate} \text{For each } y \in \Gamma(t^{n+1}), \quad \gamma(t^{n+1},y) = \frac{\gamma(t^n,(T^{-n}(y)))}{\lvert \nabla T^n (T^{-n}(y)) \tau(T^{-n}(y))\lvert},\tag{40}\] where we have used the shortcuts \[T^n = X_{t^n}(t^{n+1},\cdot) : \Gamma(t^n) \to \Gamma(t^{n+1}), \text{ and } T^{-n} = (T^n)^{-1}: \Gamma(t^{n+1}) \to \Gamma(t^{n}),\] for the flow of \(u(t,x)\) between \(t^n\) and \(t^{n+1}\) and its inverse, respectively. The update of \(\gamma(t^n,\cdot)\) into \(\gamma(t^{n+1},\cdot)\) from this formula thus requires to:

1. Calculate the position \(T^{-n}(y)\) on \(\Gamma(t^n)\) associated to each vertex \(y\) on \({\mathcal{L}}^{n+1}\);

2. Calculate the tangential derivative of the flow mapping \(T^n\) on \(\Gamma(t^n)\);

3. Assemble the expression 40 from these data.

In numerical practice, the flow mapping \(T^n\) (resp. its inverse \(T^{-n}\)) can be calculated by solving the ordinary differential equation 37 (resp. its inverse, forward expression) starting from each vertex of \(\Gamma(t^{n+1}\)) (resp. of \(\Gamma(t^n)\)). However, when evaluating 40 , it is crucial to ensure that the particle positions \(T^{-n}(y) \in \Gamma(t^n)\) attached to the vertices \(y\) of \({\mathcal{L}}^{n+1}\) belong exactly to the mesh \({\mathcal{L}}^n\) of \(\Gamma(t^n)\), since \(\gamma(t^n,\cdot)\) is solely defined on this mesh. To achieve this, we rely on the following remark about the transformation of the arc length through a diffeomorphism. Let \(\Gamma\) be a smooth 2d curve with endpoints \(c_0\), \(c_1\); the normalized arc length function \(s_\Gamma: \Gamma \to [0,1]\) reads: \[\forall x \in \Gamma, \quad s_\Gamma(x) = \frac{1}{\lvert \Gamma \lvert }\int_{\Gamma_{c_0,x}} \:\text{\rm d}s,\] where \(\lvert \Gamma \lvert\) is the length of \(\Gamma\) and we recall that \(\Gamma_{c_0,x}\) stands for the portion of curve comprised between \(c_0\) and \(x\). Let now \(T : {\mathbb{R}}^2 \to {\mathbb{R}}^2\) be a smooth diffeomorphism. This function then transforms as follows under the effect of \(T\): \[s_{T(\Gamma)}(T(x)) = \frac{1}{\lvert T(\Gamma) \lvert} \int_{T(\Gamma)_{T(c_0),T(x)}} \:\text{\rm d}s = \frac{1}{\lvert T(\Gamma) \lvert} \int_{\Gamma_{c_0,x}} \lvert \text{\rm com}(\nabla T) n \lvert \:\text{\rm d}s\] Applying this observation to the mapping \(T^{-n}\), it is possible to identify the arc length on \(\Gamma(t^n)\) associated to the particle position \(T^{-n}(y)\) on the mesh \({\mathcal{L}}^n\) of \(\Gamma(t^n)\) associated to each vertex \(y\) of \({\mathcal{L}}^{n+1}\).

Calculation of the velocity field

At each iteration \(n=0,\ldots\) of the evolution, the curve \(\Gamma(t^n)\) is discretized by a line mesh \({\mathcal{L}}^n\); once the vortex strength \(\gamma(t^n,x)\) is computed at the vertices of \({\mathcal{L}}^n\), the formula 36 for the velocity \(u(t^n,\cdot)\) is evaluated by quadrature. This task is a little delicate for two reasons:

1. The kernel \(K(x,y)\) featured in 36 blows up when \(x = y\), and this formula makes sense as a Cauchy principal value. As suggested in [64], to ease the numerical computation, we rely on the so-called “vortex-blob” approximate formula: \[\label{eq46blobpv} u(t,x) \approx \frac{1}{2\pi} \int_{\Gamma(t)} \gamma(t,y) \frac{(x-y)^\perp}{\lvert x - y \lvert^2 + \delta^2}\:\text{\rm d}s(y),\tag{41}\] where the presence of the “small” parameter \(\delta>0\) at the denominator of the above integrand removes the singularity of the kernel.

2. Often in practice, the vortex strength \(\gamma(t,x)\) takes infinite values at the endpoints of \(\Gamma(t)\) while still being an integrable function on \(\Gamma(t)\), see for instance the expression 11 used in the numerical example of [sec:sec46vortex]. To alleviate this issue, we use a quadrature formula such as the simple midpoint rule, that evaluates the integrand of 41  using quadrature points located inside the edges of \(\Gamma(t)\) and not at its endpoints.

10 Details about the treatment of the laser path optimization problem↩︎

This appendix details the calculation of the shape derivative of a slightly more general version of the functional 17 , associated to the laser path optimization example of [sec:sec46am]. We consider the quantity \[J(\Gamma) = \int_D j(x,u_\Gamma(x)) \:\text{\rm d}x,\] depending on the open curve \(\Gamma\) contained in the fixed computational domain \(D\) via the solution \(u_\Gamma \in H^1(D)\) to the boundary-value problem 16 and \(j : {\mathbb{R}}^2_x \times {\mathbb{R}}_u \to {\mathbb{R}}\) is a smooth function, satisfying suitable growth conditions. For any \(x \in D\), \(u \in {\mathbb{R}}\), we denote by \(\nabla_x j(x,u)\) and \(\frac{\partial j}{\partial u}(x,u)\) the gradient and derivative of the partial mappings \(x \mapsto j(x,u)\) and \(u \mapsto j(x,u)\), respectively.

The shape derivative of \(J(\Gamma)\) is the subject of the next proposition.

where the adjoint state \(p_\Gamma\) is the \(H^1(D)\) solution to the following boundary-value problem: \[\label{eq46pGammam} \left\{ \begin{array}{cl} -\text{\rm div}(\gamma \nabla p_\Gamma) + \beta p_\Gamma = - \frac{\partial j}{\partial u}(x,u_\Gamma)& \text{in } D, \\ \gamma \frac{\partial p_\Gamma}{\partial n} = 0 & \text{on } \partial D. \end{array} \right.\tag{42}\]

11 Mathematical treatment of the fault line detection example↩︎

This appendix provides mathematical details about the treatment of the fault line detection example of [sec:sec46fault]: the consistency of the penalized approximation of the “ideal”, discontinuous displacement field introduced in there is detailed in [sec:app46justifbroken], and the shape derivative of the least-square functional \(J(\Gamma)\) in 24 used for inversion is calculated in [sec:sec46sdfault].

Justification of the approximate model for broken spaces

In this section, we rigorously justify the approximate physical model considered in [sec:sec46fault]. For simplicity, we assume that the discontinuity line \(\Gamma\) is closed: the smooth, bounded computational domain \(D \subset {\mathbb{R}}^2\) is divided into two complementary regions \(\Omega_1, \Omega_2\), separated by the interface \(\Gamma\): \[\Omega_1 \cap \Omega_2 = \emptyset, \:\: \overline{D} = \overline{\Omega_1} \cup \overline{\Omega_2}, \quad \Gamma = \partial \Omega_1 \cap \partial \Omega_2.\] The precise situation of [sec:sec46fault], where the support \(\Gamma\) of the discontinuity in the model is an open line, is treated analogously since the imposed jump \(g_\Gamma\) on \(\Gamma\) has a smooth extension by \(0\) to a closed curve \(\widetilde{\Gamma}\) extending \(\Gamma\), i.e. \(g \in \widetilde{H}^{1/2}(\Gamma)^2\). Still for simplicity, we consider the counterpart situation of that in [sec:sec46fault] arising in the physical context of the conductivity equation: the physical behavior of the system presenting the fault is described by the solution \(u \in H^1(D \setminus \Gamma)\) to the following problem: \[\label{eq46conducfrac} \left\{ \begin{array}{cl} -\text{\rm div}(\gamma \nabla u) = 0 & \text{in } D, \\ u = 0 & \text{on } \partial D_D, \\ \gamma \frac{\partial u}{\partial n} = 0 & \text{on } \partial D \setminus \overline{\partial D_D},\\ \left[u\right] = g_\Gamma \text{ and } \left[\gamma \frac{\partial u}{\partial n}\right] = 0 & \text{on } \Gamma, \end{array} \right.\tag{49}\] where \(n\) stands for the unit normal vector to \(\Gamma\), pointing outward \(\Omega_2\), see 17 (a). Note that since \(\Gamma\) is fixed in this section, for notational simplicity, we omit the \(_\Gamma\) subscript referring to it in the solution of 49 .

According to [sec:sec46fault], the approximate version of 49 is obtained by considering two functions \(u_{1,\varepsilon}\), \(u_{2,\varepsilon}\) accounting for (approximations of) the restrictions of \(u\) to \(\Omega_1\) and \(\Omega_2\), respectively: these are both defined on \(D\) as a whole by filling the complementary regions \(D \setminus \overline{\Omega_1}\) and \(D \setminus \overline{\Omega_2}\) with a material with “soft” conductivity. Precisely, we consider the variational problem: \[\begin{gather} \label{eq46conduceappfrac} \text{Search for } (u_{1,\varepsilon}, u_{2,\varepsilon}) \in H^1_{\partial D_D}(D)^2 \text{ s.t. } \forall (v_1, v_2) \in H^1_{\partial D_D}(D)^2, \\ \int_D \gamma_{1,\varepsilon} \nabla u_{1,\varepsilon} \cdot \nabla v_1 \:\text{\rm d}x + \int_D \gamma_{2,\varepsilon} \nabla u_{2,\varepsilon} \cdot \nabla v_2 \:\text{\rm d}x + \frac{1}{\varepsilon} \int_\Gamma (u_{2,\varepsilon} - u_{1,\varepsilon}) (v_2 -v_1) \:\text{\rm d}s = \frac{1}{\varepsilon} \int_\Gamma g_\Gamma (v_2 -v_1)\:\text{\rm d}s, \end{gather}\tag{50}\] where we have introduced the conductivity coefficients: \[\gamma_{1,\varepsilon}(x) =\left\{ \begin{array}{cl} \gamma(x) & \text{if } x \in \Omega_1,\\ \varepsilon\gamma(x) & \text{if } x \in \Omega_2, \end{array} \right. \text{ and }\gamma_{2,\varepsilon}(x) =\left\{ \begin{array}{cl} \varepsilon\gamma(x) & \text{if } x \in \Omega_1,\\ \gamma(x) & \text{if } x \in \Omega_2, \end{array} \right.\]

In order to examine the consistency of the approximation 50 with the exact problem 49 as \(\varepsilon\to 0\), let us introduce the Dirichlet extensions \(u_i\) of \(u\) from \(\Omega_i\) to \(D\) for \(i=1,2\), i.e.: \[\left\{ \begin{array}{cl} u_i= u & \text{in } \Omega_i, \\ -\text{\rm div}(\gamma \nabla u_i) = 0 & \text{on } D\setminus \overline{\Omega_i}. \end{array} \right.\] Note that, as a consequence of the standard regularity theory for elliptic equations, both functions are smooth on \(\overline{\Omega_1}\) and \(\overline{\Omega_2}\): \[\label{eq46ustsmooth} \lvert\lvert u_1\lvert\lvert_{{\mathcal{C}}^{1,\alpha}(\Omega_1 \cup \Omega_2)} + \lvert\lvert u_2\lvert\lvert_{{\mathcal{C}}^{1,\alpha}(\Omega_1 \cup \Omega_2)} \:\: \leq \:\:C \lvert\lvert g_\Gamma \lvert\lvert_{{\mathcal{C}}^{1,\alpha}(\Gamma)},\tag{51}\] see e.g. §9.6 in [101]. The desired consistency result is then the following.

Calculation of the shape derivative

In this section, we detail the calculation of the shape derivative of the function \(J(\Gamma)\) at play in the shape optimization problem 24 used to reconstruct a fracture set inside a background medium.

The model context of the conductivity equation

To simplify the exposition, we first detail this derivation in the scalar context of the conductivity equation and for a generic functional \(J(\Gamma)\), of the form: \[\label{eq46Jfracgen} J(\Gamma) = \int_{D} j(x,u_\Gamma) \:\text{\rm d}x,\tag{52}\] where \(j : D \times {\mathbb{R}}\to {\mathbb{R}}\) is a given, smooth function, and \(u_\Gamma \in H^1_{\Gamma_D}(D \setminus \Gamma)\) is the solution to the conductivity equation 49 .

Our analysis starts with a lemma about the Lagrangian derivative of the normalized arc length function \(s_\Gamma : \Gamma \to [0,1]\), involved in the definition 22 of the slip function \(g_\Gamma\). Throughout this section, we denote by \(\tau:\Gamma \to {\mathbb{R}}^2\) the unit tangent vector to \(\Gamma\), oriented in such a way that for any point \(x \in \Gamma\), \((\tau(x),n(x))\) is a direct orthonormal frame of \({\mathbb{R}}^2\). In particular, \(\tau(c_0) = -n_\Sigma(c_0)\) and \(\tau(c_1) = n_\Sigma(c_1)\).

We now come to the main result of this section.

Extension in the context of the elasticity equations

We now proceed to the calculation of the shape derivative of the functional \(J(\Gamma)\) of interest, framed in the context of the linear elasticity system. Here again, to set ideas, we consider a function of the domain of the form: \[J(\Gamma) = \int_D j(x,u_\Gamma(x)) \:\text{\rm d}x,\] where \(j: D \times {\mathbb{R}}^2 \to {\mathbb{R}}\) is a smooth enough function, and \(u_\Gamma\) is the elastic displacement, solution to the system 19 . The result of interest is the following.

References↩︎