Ph J Yvonnet

D Lorong

F Ryckelynck

Chinesta

Simulating dynamic thermo-elasto-plasticity in large transformations with adaptive refinement in the natural element method : application to shear banding

Keywords: Natural element method, dynamic thermo-elastoplasticity, large transformations, adaptive refinement, shear bands

Simulation of forming processes usually involves large transformations and/or stain localization. In the context of the finite element method (FEM), a common approach to deal with such phenomena is the use of continuous adaptive remeshing. Nevertheless, constructing a quality mesh is a delicate task, especially for complex 3D problems. The natural element meshfree method is a recent numerical technique which uses the features of the Delaunay triangulation of the set of nodes, providing an accuracy equivalent to the quadrilateral/hexahedral finite elements, even if the Delaunay triangles are very distorted. In this context, inserting or removing nodes is an easy task, as no special geometrical criterion is imposed on the relative position of the nodes. Furthermore, the use of a stabilized conforming nodal integration allows to define nodal internal variables, which simplifies transfer of these data through successive updates of the reference configuration in a Lagrangian procedure. In this paper, a framework for practical implementation of the NEM in the context of explicit thermo-elastoplasticity at finite strains is provided. An adaptive strategy, based on the attractive features of the NEM is developped, including simple error indicators for history-dependent nonlinear problems, a refinement procedure for unstructured scattered of nodes based on the underlying Voronoi diagram, and the use of a nodal integration to alleviate issues associated with internal variables transfer. This promising technique for forming processes simulation is illustrated in the context of dynamic shear bands propagation after impact.

1

Introduction

Simulation of forming processes usually involves large transformations and/or strain localization. In the context of the finite element method (FEM), it is necessary to construct a mesh satisfying some quality criteria for a given deformed geometry and which also conforms the size of its elements to the localized phenomena, i.e. strain localization or cracks. Triangular and tetrahedral meshes can be automatically generated based on the Delaunay criterion in the FEM whatever the complexity of geometry shapes of the analyzed domain. Nevertheless, there are no guarantee that no flat elements can be generated, and the triangular and tetrahedral elements are rarely adopted due to their poor accuracy. Quadrilateral/hexahedral elements are usually preferred for its better accuracy, resulting in the great difficulty of mesh generation, especially for complex geometries and in the context of adaptive refinement. In recent years, meshless approximations have become interesting and promising methods in solving partial differential equations due to their flexibility in practical applications. Several meshless methods, such as, the smooth particle hydrodynamics (SPH) [START_REF] Lucy | A numerical approach to the testing of fusion process[END_REF], the diffuse element method (DEM) [START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF], the element-free Galerkin method (EFG) [START_REF] Belytschko | Element-free Galerkin methods[END_REF], the reproducing kernel particle method (RKPM) [START_REF] Liu | Reproducing Kernel Particle Methods[END_REF], the partition of unity (PUM) [START_REF] Babuska | The partition of unity finite element method: basic theory and applications[END_REF], and the natural element method (NEM) [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF] have been proposed. The NEM offers outstanding properties for the simulation of problems involving large transformations and adaptive refinements due to the following properties : (a) essential boundary can be enforced directly as the NEM shape functions satisfies the Kronecker delta property, and are also strictly linear over the boundaries of the convex hull (it has been extended to non convex geometries in [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM)[END_REF]); (b) the support of the shape functions, constructed on the basis of the underlying Voronoi diagram (dual of the Delaunay triangulation), automatically adapts its size and shape to the local neighborhood, whatever the complexity and density of of the surrounding nodal distribution; (c) the accuracy does not depend significantly on the relative position of the nodes (shape of the Delaunay triangles), and it reslults equivalent to quadrangular/hexahedral finite elements [START_REF] Yongchang | A meshless local natural neighbor interpolation method for stress analysis of solids[END_REF] even when Delaunay triangles become very distorted. Is noteworthy that this property greatly simplifies the treatment of complex three-dimensional domains, as the Delaunay triangulation is unique for a given set of nodes.

In this paper, we analyze the application of the C-NEM to the simulation of dynamic thermo-elastoplastic behaviors, involving large transformations combined with an adaptive strategy based on the attractive features of the NEM. In section 2, the natural element method is briefly reviewed, as well as its extension for arbitrary (non convex) geometries that results in the constrained natural element method (C-NEM) [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM)[END_REF]. In section 3, a Lagrangian formulation of the coupled thermo-mechanical problem in large transformations is established. In section 4, a strategy for refinement and adaptivity in the NEM (or C-NEM) context is proposed, based on the attractive features of the NEM. A simple error indicator using the Zienkiewicz-Zhu [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF] paradigm in tandem with the stabilized conforming nodal integration is proposed. A refinement strategy based on the underlying Voronoi diagram is proposed, which allows refinement through arbitrary unstructured clouds of nodes. A nodal-NEM formulation [START_REF] Yoo | Stabilized conforming nodal integration in the naturalelement method[END_REF] is also used to define the internal variables, usually defined at the integration points in the context of FEM, at the nodes, which greatly simplifies the transfer of these variables through successive updates of the reference configuration. Finally, the technique is applied in section 5 in the context of dynamic shear bands propagation in metals.

The constrained natural element method

Natural neighbor interpolation

We briefly touch upon the foundation of Sibson's natural neighbor coordinates (shape functions) that are used in the natural element method. For a more in-depth discussion on the Sibson interpolant and its application for solving second-order partial differential equations, the interested reader can refer to Braun and Sambridge [START_REF] Sambridge | Geophisical parameterization and interpolation of irregular data using natural neighbors[END_REF], and Sukumar et al. [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]. The NEM interpolant is constructed on the basis of the Voronoi diagram. The Delaunay tessellation is the topological dual of the Voronoi diagram.

Consider a set of nodes S = {n 1 , n 2 , . . . , n N } in dim . The Voronoi diagram is the subdivision of dim into regions T i (Voronoi cells) defined by:

T i = {x ∈ dim : d(x, x i ) < d(x, x j ), ∀j = i}, ∀ i (1)
The Sibson coordinates of x with respect to a natural neighbor n i (see Fig. 1) is defined as the ratio of the overlap area (volume in 3D) of their Voronoi cells to the total area (volume in 3D) of the Voronoi cell related to point x:

φ i (x) = Area(af ghe) Area(abcde) (2)
If the point x coincides with the node n i , i.e. (x = x i ), φ i (x i ) = 1, and all other shape functions are zero, i.e. φ i (x j ) = δ ij (δ ij being the Kronecker delta). The properties of positivity, interpolation, and partition of unity are then verified [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]:

⎧ ⎨ ⎩ 0 ≤ φ i (x) ≤ 1 φ i (x j ) = δ ij n i=1 φ i (x) = 1 (3)
The natural neighbor shape functions also satisfy the local coordinate property [START_REF] Sibson | A vector Identity for the Dirichlet tessellations[END_REF], namely:

x = n i=1 φ i (x)x i (4)
which combined with Eq. ( 3) implies that the natural neighbor interpolant spans the space of linear polynomials (linear completeness).

It turns out that the support of φ i (x) is the union of the n circles (spheres in 3D) passing through the vertices of the n Delaunay triangles (tetrahedra) containing the node n i (in this case n is the number of natural neighbors of node n i ). The support of a node n i in a particular nodal distribution is depicted in figure 2.

Natural neighbor shape functions are C ∞ at any point except at the nodes, where they are only C 0 , and on the boundary of the Delaunay circles (spheres in 3D) where they are only C 1 , because of the discontinuity in the neighbors nodes across these boundaries. Hiyoshi and Sugihara [START_REF] Hiyoshi | Improving continuity of Voronoi-based interpolation over Delaunay spheres[END_REF] have shown that the Sibson interpolant belongs to a more general class of Voronoi-based interpolants, called k-th order standard coordinates, proving that the interpolant generated by the k-th order standard coordinates have C k continuity on the Delaunay circles (spheres) boundaries. In this context, the Sibsonian and non-Sibsonian (Laplace) coordinates [START_REF] Belikov | The non-Sibsonian interpolation : a new method of interpolation of the values of a function on an arbitrary set of points[END_REF] results to be the standard coordinates of order 1 and 0, respectively.

Another important property of this interpolant is the ability to reproduce linear functions over the boundary of convex domains. The proof can be found in Sukumar et al. [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]. An illustration is depicted in Fig. 1 (b): as the areas associated to points on the boundary become infinite, the contribution of internal points vanish in the limit when the point approaches the convex boundary, and the shape functions associated with nodes n 1 and n 2 become linear on the segment (n 1n 2 ). This is not true in the case of non convex boundaries, and the next section focuses on an approach to circumvent this difficulty.

Consider an interpolation scheme for a function (vectorial, scalar or tensorial) u(x), in the form: where u i are the nodal values at the n natural neighbor nodes, and φ i (x) are the shape functions associated with each neighbor node. It is noted that Eq. ( 5) defines a local interpolation scheme. Thus, the trial and test functions used in the discretization of the variational formulation describing both the mechanical and thermal problems treated in this paper take the form of Eq. ( 5).

u h (x) = n i=1 φ i (x) u i (5)

The Constrained natural element method

Constrained Voronoi diagram

In its original form [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF], the NEM can only be applied to strictly convex domains. For non-convex domains, two main issues occur : (a) It was proved in [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF][START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM)[END_REF] and [START_REF] Cueto | Imposing essential boundary conditions in the natural elements method by means of density-scaled alpha-shapes[END_REF] that a loss of linearity in the interpolation along boundaries of non convex domains appear. Thus essential boundary conditions can only be imposed directly over convex boundaries; (b) For strongly non-convex domains (cracks, auto-contact...) some spurious influences between nodes of the boundaries appear [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM)[END_REF]. An additional treatment is thus required to maintain all the properties of the NEM for any geometry.

In order to avoid these drawbacks, we have proposed in a previous paper [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM)[END_REF] an extension of the NEM in which a visibility criterion is introduced in order to restrict influent nodes among natural neighbors. The computation of the shape functions is done on the basis of the so-called constrained (or extended) Voronoi diagram (CVD), which is the strict dual to the constrained Delaunay triangulation, introduced by Seidel in [START_REF] Seidel | Constrained Delaunay triangulations and Voronoi diagrams with obstacles In "1978-1988[END_REF], instead of the Voronoi diagram (see [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM)[END_REF] for further details). The intersection between the CVD and the domain results into new cells T C i , called constrained Voronoi cells, defined formally by:

T C i = {x ∈ n : d(x, x i ) < d(x, x j ), ∀j = i, S x→n i ∩ Γ = ∅, S x→n j ∩ Γ = ∅} (6)
where Γ is the domain boundary, composed by a set of segments l i ∈ L and S a→b denotes the segment between the points a and b. In this framework, a point located inside a cell T C i is closer to the node n i than to any other visible node n j .

The constrained Delaunay triangulation does not always exist in 3D without adding new nodes, as shown in [START_REF] Schönhardt | Uber die zerlegung von dreieckspolyedern in tetraeder[END_REF]. Nevertheless, some techniques for constructing 3D constrained Delaunay tessellations are available and provided in [START_REF] Shewchuck | Tetrahedral mesh generation by delaunay refinement[END_REF][START_REF] Shewchuck | Sweep algorithms for constructing higher-dimensional constrained Delaunay triangulations[END_REF] by adding Steiner points.

The constrained natural element approximation

In order to solve partial differential equations defined on non convex domains, or to reproduce functional discontinuities, we consider the following approximation for both the trial and the test functions:

u h (x) = V i=1 φ C i (x)u i ( 7
)
where V is the number of natural neighbors visible from point x and φ C i is the constrained natural neighbor shape function related to the i-th node at point x. The computation of the C-n-n (constrained natural neighbor) shape functions is similar to the natural neighbor shape function, when one proceed using the constrained Voronoi diagram introduced previously. It was shown in [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM)[END_REF] and [START_REF] Yvonnet | The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces[END_REF] that the use of the constrained Voronoi diagram does not affect the properties of the NEM interpolation, allowing the extension of the linearity of the shape functions to any geometry, convex or not.

The ability of the C-NEM for treating problems involving cracks has been illustrated in [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM)[END_REF] and in the context of moving interfaces in [START_REF] Yvonnet | The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces[END_REF], where we have shown how the C-NEM simplifies the treatment of material discontinuities in meshfree methods, due to the continuity of the approximation across interfaces (consequence from the interpolant character and linearity of the C-NEM shape functions over any external boundary or internal interface). It is noteworthy that no size parameter is involved in the definition of the influent nodes (neighbors), which is essential for the robustness of refinement procedures, in which the nodal density varies in some parts of the domain.

Formulation of the coupled thermo-mechanical problem

Preliminaries

In this section, we summarize briefly some fundamental kinematic relations and introduce the appropriate notations. We consider two configurations of a body B : the first one, the reference configuration Ω 0 ⊂ dim with external boundary ∂Ω 0 (not necessarily the initial configuration), where X denotes the coordinates of a point in this configuration. The second, called current configuration, is denoted by Ω x ⊂ dim with external boundary ∂Ω x at time t, with x the coordinates of a point in the current configuration. x is related to X by :

x = X + u(X, t) (8)
the deformation gradient is defined by:

F = ∂x ∂X = 1 + ∇ X u, J = det(F) > 0 (9)
From now on, we use a hyperelastic formulation based on the multiplicative decomposition of the deformation gradient, with elastic response described by a hyperelastic stored energy function. In this context, elastic predictor becomes exact and the need for incrementally objective algorithms is entirely avoided.

The basic hypothesis underlying this approach to finite strain elastoplasticity is the multiplicative split of the deformation gradient, F, into elastic and plastic parts:

F = F e F p (10)
This assumption, firstly proposed by Lee [START_REF] Lee | Elastic-plastic deformation at finite strains[END_REF], admits the existence of a local unstressed intermediate configuration. Following the multiplicative split of F, the velocity gradient, L = ḞF -1 can be decomposed additively as

L = L e + L p (11)
where L e and L p are, respectively, the elastic and plastic contributions defined by:

L e = Ḟe [F e ] -1 , L p = F e Ḟp [F p ] -1 [F e ] -1 (12)
Similarly, the stretching tensor, D := sym[L], can be decomposed as:

D = D e + D p (13)
with the elastic and plastic stretching tensors given by

D e = sym[L e ], D p = sym[L p ] (14)

Hyperelastic constitutive law

From polar decomposition, F e is given by

F e = R e U e (15)
where U e and R e are, respectively, the elastic right stretch tensor and the elastic rotation.

Let e denote the Eulerian (or spatial) logarithmic strain tensor [3]

e = ln[U e ] (16)
where ln[.] above denotes the tensor logarithm of (.) which involves spectral decomposition of U e . Following Peric et al. [START_REF] Perić | On adaptive strategies for large deformations of elastoplastic solids at finite strains: computational issues and industrial applications[END_REF], we assume the existence of a quadratic strain energy function ψ e ( e ) in the form of a scalar symmetric function of its stretches λ i (i = 1, 2, 3) given by:

ψ e (λ e 1 , λ e 2 , λ e 3 ) = μ ln(λ e 1 ) 2 + ln(λ e 2 ) 2 + ln(λ e 3 ) 2 + 1 2 λ(J e ) 2 (17)
where μ and λ are Lamé's parameters and (J e ) = λ e 1 λ e 2 λ e 3 is the Jacobian. After applying standard procedure, the following hyperelastic constitutive equation is obtained:

T = ∂ψ e ∂ e = C e : e ( 18
)
where T is the rotated stress tensor. Assuming incompressibility of the plastic flow, it is expressed as:

T = [R e ] T τ R e (19)
where τ = Jσ is the Kirchhoff stress tensor, C e is the fourth-order isotropic elastic tensor. Further details about the thermo mechanical foundations can be found in [START_REF] Eterovic | A hyperelastic-based large strain elasto-plasstic constitutive formulation with combined isotropic-kiematic hardening using the logarithmic stress and strain measures[END_REF].

Numerical integration of the constitutive equations

Typically, within an incremental numerical procedure for solving history dependent problems, a numerical approximation to the material constitutive law is needed to update stresses τ as well as the internal variables α ≡ ¯ p (equivalent plastic strain) within each time (load) increment. In the present context, given the values of the variables {τ n , F p n , α n } at the beginning of a generic increment [t n , t n+1 ], an algorithm is required to update τ n+1 , F p n+1 , α n+1 at the end of the increment. Under the assumption of elastic isotropy, the elastic Eulerian logarithmic strains tensor is updated according to:

e n+1 = e trial n+1 -ΔγN n+1 (20)
which has the same format that the standard return mapping used in the infinitesimal theory [START_REF] Eterovic | A hyperelastic-based large strain elasto-plasstic constitutive formulation with combined isotropic-kiematic hardening using the logarithmic stress and strain measures[END_REF][START_REF] Simo | Algorithms for static and dynamic multiplicative plasticity that preserve the classical return-mapping algorithm schemes of the infinitesimal theory[END_REF][START_REF] Cuitiño | A material-independent method for extending stress-update algorithms from small-strain plasticity to finite plasticity with multiplicative kinematics[END_REF][START_REF] Simo | Associative coupled thermoplasticity at finite strains : formulation numerical analysis and implementation[END_REF].

N n+1 = 3 2 dev(T n+1 ) J 2 (T n+1 ) (21)
where J 2 (.) is the second invariant of (.). The trial elastic logarithmic strain e trial n+1 is given by:

e trial n+1 = ln[U e trial n+1 ] (22)
where U e trial n+1 results from polar decomposition of F e trial n+1 defined by:

F e trial n+1 = F n+1 (F p n ) -1 (23)
From Eq. ( 18) we have:

dev(T n+1 ) = 2μdev( e n+1 ) (24)
Using Eq. ( 20), Eq. ( 21) and Eq. ( 24) we obtain:

dev(T n+1 ) = dev(T trial n+1 ) -3μΔγ dev(T n+1 ) J 2 (T n+1 ) (25)
with:

T trial n+1 = C e : e trial n+1 (26)
which leads to :

J 2 (T n+1 ) = J 2 (T trial n+1 ) -3μΔγ (27)
We assume that f y trial n+1 = J 2 (T trial n+1 )σ y (Δγ) is a non-linear scalar function. In the case of a plastic increment (f y trial n+1 > 0), we must solve for Δγ :

J 2 (T trial n+1 ) -3μΔγ -σ y (Δγ) = 0 ↔ J 2 (T n+1 ) -σ y (Δγ) = 0 (28)
A classical Newton-Raphson procedure has been used in this work to solve the above equation.

Using Eq. ( 25) we obtain :

dev(T n+1 ) = dev(T trial n+1 ) 1 + 3μΔγ J 2 (T n+1 ) (29)
with J 2 (T n+1 ) defined accorind to Eq. ( 27) and:

T n+1 = dev(T n+1 ) + 1 3 T r(T trial n+1 ) (30)
The Cauchy stress tensor may thus be obtained as:

σ n+1 = 1 det(F n+1 ) R etrial n+1 -T T n+1 R etrial n+1 -1 (31)
Finally using Eq. ( 20) and Eq. ( 22), and assuming that :

R e n+1 = R e trial n+1 (32)
Box 1. Algorithm for integration of constitutive equations. (i) For given displacement u n+1 , evaluate total deformation gradient

F n+1 = 1 + ∇ X u n+1 = 1 + B(X)u n+1
with B(X) a matrix containing the shape functions derivatives in the reference configuration. (ii) Evaluate elastic trial deformation gradient

(F e ) trial n+1 = (F n+1 )(F p n ) -1 (
J 2 (T n+1 ) -σ n+1 y (Δγ) = 0, with J 2 (T n+1 ) = J 2 (T trial n+1 ) -3μΔγ (viii) Update Cauchy stress tensor σ n+1 = J -1 n+1 [R e n+1 ] -T T n+1 R e -1 n+1 , with T n+1 = dev(T n+1 ) + 1 3 T r(T trial n+1 )1 dev[T n+1 ] = dev[T trial n+1 ] 1+Ξ , Ξ = 3μΔγ J 2 (T n+1 ) (ix) Update plastic part of deformation gradient F p n+1 = exp[ΔγN n+1 ]F p n , N n+1 = 3 2 dev(T n+1 ) J 2 (T n+1 )
we obtain, after some calculations, the following incremental law for F p (see [START_REF] Eterovic | A hyperelastic-based large strain elasto-plasstic constitutive formulation with combined isotropic-kiematic hardening using the logarithmic stress and strain measures[END_REF][START_REF] Simo | Algorithms for static and dynamic multiplicative plasticity that preserve the classical return-mapping algorithm schemes of the infinitesimal theory[END_REF][START_REF] Cuitiño | A material-independent method for extending stress-update algorithms from small-strain plasticity to finite plasticity with multiplicative kinematics[END_REF]) :

F p n+1 = exp[ΔγN n+1 ]F p n ( 33
)
As noticed in [START_REF] Perić | On adaptive strategies for large deformations of elastoplastic solids at finite strains: computational issues and industrial applications[END_REF], as a consequence of the exponential mapping in the implicit integration of the plastic flow rule, the incompressibility of the plastic flow for pressure insensitive yield criteria is carried over exactly to the incremental rule [START_REF] Shewchuck | Tetrahedral mesh generation by delaunay refinement[END_REF]. The algorithm, therefore, generalizes the standard return mapping algorithms [START_REF] Simo | Computational inelasticity[END_REF] of the infinitesimal theory. The overall algorithm for the incremental stress update is outlined in Box 1.

Explicit Lagrangian procedure

With the principle of virtual work as a basis of kinematically based C-NEM solution scheme, the corresponding continuum incremental boundary value problem is formulated in the spatial configuration as follows.

Ω t ρ(t)ü•ηdΩ t + Ω t σ t : ∇ x ηdΩ t = Ω t ρ(t)b•ηdΩ t + ∂Ω t σ t•ηdΓ t ∀η ∈ ϑ (34)
where ϑ is the space of virtual displacements. The properties dΩ t = J t dΩ 0 and ρ 0 dΩ 0 = ρ(t)dΩ t are used, which leads to:

Ω 0 ρ 0 ü•ηdΩ 0 + Ω 0 P t : ∇ X ηdΩ 0 = Ω 0 ρ 0 b•ηdΩ 0 + ∂Ω t σ t•ηdΓ t ∀η ∈ ϑ (35)
where P denotes the first Piola-Kirschhoff stress tensor related to σ by P = JF -1 σ.

The C-NEM discretization [START_REF] Camacho | Computational modeling of impact damage in brittle materials[END_REF] of the variational form [START_REF] Seidel | Constrained Delaunay triangulations and Voronoi diagrams with obstacles In "1978-1988[END_REF] results in the discrete set of algebraic time dependent equations which may be expressed, in matrix form, as:

Mü n+1 (t) = F ext n (t) -F int n (u n , t) ( 36
)
where t is the time, M denotes the mass matrix, F int n (u, t) the internal force vector, while F ext n (t) is the external force vector, expressed, respectively, by:

M = Ω 0 ρ 0 φ T (X)φ(X)dΩ 0 (37)
F int n (u n , t) = Ω 0 J n F -1 n σ n B(X)dΩ 0 ( 38
)
F ext n = ∂Ω t σ φ T (x)tdΓ t (39)
with φ(X) a matrix containing the shape functions in the reference configuration and B(X) a matrix containing the shape functions derivatives also in the reference configuration. As shown in the next section, the use of the SCNI quadrature [START_REF] Chen | Non-linear version of stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] results in a M matrix diagonal, whose diagonal terms are given by m i = ρ 0 Ω i , with Ω i the area (volume in 3D) of the Voronoi cell related to node n i .

The velocity v = u and acceleration ü = v are approximated by using central differences with variable time steps. Thus, we have:

v n+1/2 = v n-1/2 + Δt 1 + Δt 2 2 ün (40)
u n+1 = u n + Δt 2 v n+1/2 (41)
Finally, the general explicit algorithm is outlined as follows. Being known the initial conditions or the computed solution at time: t n : u n , v n-1/2 , ün , F n , F p n :

(i) Update displacements and velocity

v n+1/2 = v n-1/2 + Δt 1 +Δt 2 2 ün u n+1 = u n + Δt 2 v n+1/2 (ii) Evaluate σ n+1 , F p n+1 by using Box 1. (iii) Update accelerations ün+1 ( üi ) n+1 = 1 m i (F ext i ) n+1 -F int i n+1
Remark: In case of contact, ün+1 is only used as a predictor phase, which has to be corrected according to the prescribed displacements or traction.

Thermo-mechanical coupling

The weak form of the heat balance can be expressed as

Ω t ρ(t)c(t) Ṫ ηdΩ t + Ω t k(t)∇ x T • ∇ x ηdΩ t = Ω t r t ηdΩ t + ∂ 2 Ω qηdΓ t ∀η ∈ V T (42)
where c(t) is the specific heat, k(t) is the thermal conductivity for isotropic conduction, and r t a heat source related to the inelastic deformations, defined at time t. In the following, we assume c and k constant in time. V T is the space of virtual temperatures and q represents the heat transfer at the boundary ∂ 2 Ω. Following similar arguments as in the mechanical problem [START_REF] Shewchuck | Sweep algorithms for constructing higher-dimensional constrained Delaunay triangulations[END_REF], the weak form of the heat balance can be expressed in the reference configuration by

Ω 0 ρ 0 c Ṫ ηdΩ 0 + Ω 0 k∇ X TF -T F -1 ∇ X ηdΩ 0 = Ω 0 det(F)r t ηdΩ 0 + + ∂ 2 Ω qηdΓ t ∀η ∈ V T (43)
The C-NEM discretization of Eq. ( 43) results in the ODE system:

C Ṫ + KT = Q (44)
which in the context of an explicit scheme can be written as

CT n+1 = [C + ΔtK]T n + ΔtQ n with C = Ω 0 ρ 0 cφ T (X)φ(X)dΩ 0 (45)
K = Ω 0 kB T (X)F -T F -1 B(X)dΩ 0 (46)
Q n = ∂Ω t 2 φ T (x)qdΓ t (47)
The matrix C becomes diagonal in the context of C-NEM when stabilized conforming nodal integration is used, being the diagonal terms c i = ρ 0 cΩ i . The heat source resultig from the inelastic deformations is given by

r t = χσ t : Dp t (48) with Dp t = sym Ḟp t (F p t ) -1 (49)
where σ t is the Cauchy stress tensor et time t, and χ is the Taylor-Quinney parameter [START_REF] Taylor | The latent energy remaining in a metal after cold working[END_REF] representing the fraction of plastic work converted into heat. In the present work, we have used χ = 0.9 [START_REF] Marusich | Modelling and simulation of high-speed machining[END_REF]. The thermo-mechanical coupling is carried out by the effects of inelastic deformations [START_REF] Zienkiewicz | Recovery procedures in error estimation and adaptivity Part I: Adaptivity in linear problems[END_REF] and by the softening effects due to temperature in the hardening law. The hardening law adopted in the present simulation is given by the classical Johnson-Cook law [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF]

σ y (¯ p ) = [A + B(¯ p ) n ] 1 + Cln ˙ p ˙ p 0 1 -T -T 0 T f -T 0 m (50)
where A, B, C, m, n and ˙ p 0 are material parameters, ¯ p represents the equivalent plastic strain and ˙ p represents the rate of plastic strain. T and T 0 correspond to the current and initial temperature respectively.

When mechanical and thermal problems are coupled, a staggered solution approach is usually adopted, which solves the mechanical and the thermal problems in an uncoupled manner with data exchange performed at the end of each time step or increment. In particular, the nodal temperatures are transferred to the mechanical procedure, while plastic work is communicated to the thermal solver.

Refinement and adaptivity in the C-NEM

In this section, an adaptive strategy is proposed for non-linear problems with history-dependent internal variables in the C-NEM context, including : (a) a Zienkiewicz-Zhu error indicator [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF] based on equivalent plastic strain; (b) a transfer technique based on the stabilized conforming nodal integration [START_REF] Chen | Non-linear version of stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] and (c) a refinement strategy for domains evolving in time based on the Voronoi cells.

Transfer of internal variables with nodal integration

In the context of the finite element method, the transfer of internal variables between successive remeshing operations is an important issue (see i.e. [START_REF] Bathe | Error indicators and adaptive remeshing in large deformation analysis[END_REF][START_REF] Camacho | Computational modeling of impact damage in brittle materials[END_REF][START_REF] Perić | On adaptive strategies for large deformations of elastoplastic solids at finite strains: computational issues and industrial applications[END_REF]). In this work, we use the stabilized conforming nodal integration proposed by Chen et al. [START_REF] Chen | Non-linear version of stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] to define all variables at the nodes, in order to avoid projection between two successive actualization of the reference configuration (update of the Voronoi diagram and shape functions).

We thus consider a deformation gradient at node n i :

Fi = 1 |Ω i | Ω i ∇ X u h (x) + 1 dΩ = 1 + ∇X u h (x i ) = 1 + Bi u i (51)
where Ω i is a representative domain around the node (typically a Voronoi cell) (see fig. 3). Introducing C-NEM discretization scheme into (51) we obtain:

∇X u h (x i ) = ⎡ ⎢ ⎢ ⎣ u h ,X (x i ) u h ,Y (x i ) v h ,X (x i ) v h ,Y (x i ) ⎤ ⎥ ⎥ ⎦ = Bi u i ( 52
)
where Bi is expressed by:

Bi = ⎡ ⎢ ⎢ ⎣ φ1,X 0 φ2,X 0 . . . φN,X 0 φ1,Y 0 φ2,Y 0 . . . φN,Y 0 0 φ1,X 0 φ2,X . . . 0 φN,X 0 φ1,Y 0 φ2,Y . . . 0 φN,Y ⎤ ⎥ ⎥ ⎦ (53) with φi,X = Ω φ i,X (X)dΩ (54) φi,Y = Ω φ i,Y (X)dΩ (55)
where u i are the nodal displacements. In the context of a Lagrangian procedure, different options can be considered: 1) A total Lagrangian procedure, where the Voronoi diagram as well as the shape functions are computed only once at the beginning of the simulation;

2) An updated Lagrangian procedure where the Voronoi diagram and the shape functions are adapted every time steps.

3) A periodically updated Lagrangian procedure where the Voronoi diagram and the shape functions are computed after several total Lagrangian steps.

In the context of a total Lagrangian procedure, no internal variable transfer is necessary. Nevertheless, in updated Lagrangian procedures (2) and (3), it is necessary to update the reference configuration, which implies to reconstruct the new Voronoi diagram associated with the new nodal distribution, as well as the new shape functions. In that case, the use of the assumed gradient defined in (51) allows the definition of all internal variables at the nodes. In the case of an explicit procedure, all quantities defined in the former time step are defined at the nodes in the same way. If the same cloud of nodes is used through the whole simulation, the transfer is thus direct as integration points coincide with the nodes. If a different cloud of node is used during the simulation (due to the refinements or nodal repositioning, a direct C-NEM interpolation [START_REF] Camacho | Computational modeling of impact damage in brittle materials[END_REF] can be performed to define the internal variable the new nodes.

Error indicator based on equivalent plastic strain

In the context of the finite element method, adaptive strategies have been extensively applied and continuously developed for linear [START_REF] Ladevèze | Non Linear Computational Structural Mechanics Springer[END_REF] [1] [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF] and some class of nonlinear problems and history-dependent nonlinear problems over the last two decades or so (see i.e. Ladevèze et al. [START_REF] Ladavèze | Accuracy of elastoplastic and dynamic analysis[END_REF], Belytschko et al. [START_REF] Belytschko | Fission-fusion adaptivity in finite elements for nonlinear dynamics of shells[END_REF], Ortiz et Quigley [START_REF] Ortiz | Adaptive mesh refinement in strain localization problems[END_REF], Gallimard et al. [START_REF] Gallimard | Error estimation and adaptivity in elasto-plasticity[END_REF], among many others).

The present section aims to propose a simple error indicator for adaptive solutions of large elasto-plastic transformations in the C-NEM framework. In the context of the finite element method, an elementary procedure for the error estimation may be defined by substituting the exact solution by some post-processed field obtained from the available FEM solution. When the finite element solution is accurate enough, the post-processed solution is expected to be more accurate than the original FEM solution. In particular, the a posteriori error estimation procedure originally proposed and used by Zienkiewicz and Zhu [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF] for linear elliptic problems is based on the observation that exact stress σ may be represented accurately by smoothed stress σ * obtained by a suitable projection of stresses σ h .

In the context of meshless methods, the field σ h is smooth in the general case, unlike in the FEM. Many options have been considered to construct error indicators in meshfree methods, mostly using some recovery Zienkiewicz-Zhu fields (see the work by Liu et al. [START_REF] Liu | Enrichment of the finite element method with the reproducing kernel particle method[END_REF], You et al. [START_REF] You | Filters, Reproducing Kernel, and Adaptive Meshfree Method[END_REF], Chung et Belytschko [START_REF] Chung | An error estimate in the EFG method[END_REF], Lee et Zhou [START_REF] Lee | On error estimation and adaptive refinement for element free galerkin method: Part II:adaptive refinement[END_REF][START_REF] Lee | On error estimation and adaptive refinement for element free galerkin method: Part I: stress recovery and a posteriori error estimation[END_REF] and Lu et Chen [START_REF] Lu | Adaptive meshfree particle method[END_REF] for an overview of the recent proposed techniques). In most of the referred papers in which a MLS approximation is used, an important issue is the influence of the shape function support size on the efficiency of error estimates [START_REF] Chung | An error estimate in the EFG method[END_REF]. Furthermore, it has been shown in [START_REF] Lee | On error estimation and adaptive refinement for element free galerkin method: Part I: stress recovery and a posteriori error estimation[END_REF] that adaptive refinement in MLS meshless techniques is a delicate task due to the necessity to adapt locally the shape function support according to the local nodal density. In the following, an error indicator based on the NEM shape functions is proposed to circumvent the just referred difficulties, as the NEM shape functions support automatically adapts its shape and size to the surrounding neighborhood.

We propose a simple error indicator using the Zienckiewicz-Zhu idea in tandem with the stabilized conforming nodal integration. Let α(X, t) a variable either associated with the spatial derivatives (i.e. F(X, t)), or with internal variables (i.e. ¯ p (X, t)). In the C-NEM context, the stabilized conforming nodal integration scheme proposed by Chen et al. [START_REF] Yoo | Stabilized conforming nodal integration in the naturalelement method[END_REF] produces constant piece-wise fields associated with the different derivatives, discontinuous across the Voronoi cells. This is a consequence of the stabilization scheme used in (51). In addition, the different internal variables can be considered at the nodes, as we have nodal integration. Constant fields α h (X i , t) associated with variables α(X, t) can thus be considered in each Voronoi cell Ω i . A simple solution for recovery fields α * (X, t) with assumed better accuracy is to interpolate the nodal values of α h (X i , t) with the C-NEM interpolation scheme:

α(X, t) * = V i=1 φ C i (X)α h (X i , t) (56)
Then, the error based on the equivalent plastic strain for any cell Ω i may be defined in the reference configuration, as:

|e| 2 Ω i = Ω i0 ([¯ p (X, t)] * -[¯ p (X i , t)] h ) 2 J i dΩ 0 (57)
with:

[¯ p (X, t)] * = V j=1 φ c j (X)[¯ p (X j , t)] h , J i = det( Fi ) (58)
Where [¯ p (X j , t)] h is the equivalent plastic strain associated with the natural neighbor n j of point x, assumed constant in the constrained Voronoi cell Ω j . In order to evaluate (57), the constrained Voronoi cells are triangulated and standard Gauss quadrature is applied on the triangles.

Any kind of error indicators can be constructed on the above framework, i.e. error indicators based on rate of plastic work or based on damage [START_REF] Perić | On adaptive strategies for large deformations of elastoplastic solids at finite strains: computational issues and industrial applications[END_REF]. For the sake of simplicity, we only focus in this study on the error indicator described below. The global error is obtained, in a standard way, by :

|e| 2 Ω = N i |e| 2 Ω i .
In addition, the relative error η is defined as

η i := |e| Ω i /|e| Ω

Refinement strategy based on the Voronoi cells

To adapt the nodal density to the evolution of plastic deformation, we must define a limit global error η and a local (in each cell) limit error ηi satisfying:

η2 = N i=1 (η i ) 2 (59)
If we assure a uniform limit error ηi = η * , ∀i then:

η * = η √ N (60)
Following Zienckiewicz et. al [START_REF] Zienkiewicz | Recovery procedures in error estimation and adaptivity Part I: Adaptivity in linear problems[END_REF], we propose to determine the characteristic length of the new cell compared to its parent cell (which has been divided to generate new cells). For this purpose, we assume a rate convergence in energy norm of p (O((h p ))). We thus have:

h new i h old i = η * η i 1 p = η|e| Ω √ N |e| Ω i 1 p (61)
As noticed in former studies [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM)[END_REF][START_REF] Sukumar | The natural elements method in solid mechanics[END_REF], the convergence rate in the natural element method is of the same order as that obtained in the FEM using linear finite elements (due to the linear consistency of the Sibson interpolant), that is p = 1. To determine when a cell has to be subdivided, we assume that the characteristic length of the offspring cells are twice lower (see fig. ( 4)), for sake of simplicity. The criterion for Voronoi cell subdivision is then defined by:

h new h old ≤ 1 2 (62)
η|e| Ω √ N |e| Ω i 1 p ≤ 1 2 (63)
Figure 4 shows different refinement possibilities based on the Voronoi diagram. On the left is depicted the Voronoi diagram before insertion of new nodes. On the right is depicted the updated Voronoi diagram. In the strategy depicted in figure 4 (a), new nodes are added on the vertex of the Voronoi cell when the criterion (63) is satisfied. This strategy has been used by several authors [START_REF] You | Filters, Reproducing Kernel, and Adaptive Meshfree Method[END_REF] [START_REF] Lu | Adaptive meshfree particle method[END_REF] for refinement in unstructured scattered of nodes. In the strategy depicted in figure 4 (b), new nodes are added between the node associated with the cell and the neighbors, when the criterion (63) is satisfied. We recommend to use the strategy (b), as in the case of strategy (a), for regular grids, several additional nodes can be inserted at the same position. These two strategies allow refinement in unstructured set of nodes and are consistent with the error indicators based on the Voronoi cells described below. In addition, no parameter is associated with the shape functions support size (see Fig. 2), which simplifies the adaptive analysis when the nodal density is significantly different in some parts of the domain.

Dynamic refinement strategy

We propose the following strategy for dynamic refinement: starting for the reference configuration (at time t 0 ), several Lagrangian increments are performed until a chosen time t a . At time t a , errors are evaluated using (57). If maximal accepted error is reached in some cells, new nodes are added according to the refinement strategy described in fig 4 (b) in the reference configuration. The computation then restart from time t 0 with the new nodes until time t a . The procedure continues while the error is lower than the prescribed limit in all the cells of the domain. Then the configuration is updated at time t a , which will constitute the new reference configuration.

Numerical example

High-speed shearing is a new metal cutting operation which offers significant advantages compared to traditional shearing. In this process, the high velocity of the punch induces adiabatic shear bands which combined with the damage phenomena produces high-quality of the cutting zone, without any burr formation. This process allows to perform shearing in high strength material such as titanium or stainless steels, or to produce holes in plates with relatively large thickness (up to 20mm) in traditional metals such as aluminum.

The example shows a first insight into the simulation of a plane-strain, high-speed shearing operation using the error indicator and the refinement procedure described previously. The geometry of the problem is depicted in Fig. 5. The following material data for stainless steel 304L have been used: A = 253 MPa, B = 685 MPa, C = 0.0973 n = 0.312, m = 2.044, ˙ p 0 = 1; T 0 = 296 K, T f = 1698 K, E = 210 GPa, ν = 0.33, c = 500 J/kgK, k =5.86 J/JgmK, ρ = 7850 kg/m 3 . The initial velocity of the punch is taken as V = 11 m.s -1 .

The main aim was to illustrate the capability of the refinement procedure to capture the details of the shear bands when different geometries of the tools are used. The high speed of the punch, combined to the material characteristics, induces high strain localization. The high rate of energy generation induced by the plastic deformation and the low thermal diffusivity of the material causes localized heating, which in turn, causes thermal softening and consequently localized plastic deformation. Two test have been performed using the following geometrical parameters: R p = 5 mm, L = 20 mm, H = 10 mm and j = 0.2 mm. In a first test, zero corner radii R 1 and R 2 have been used. Figure 6 shows the refined nodal distribution during the shear band formation. Fig. 7 illustrates the evolution of the underlying Voronoi diagram used for the refinement procedure in half of the workpiece. Fig. 8 depicts the equivalent plastic strain. In this first example, the direction of the shear band is observed parallel to the punch speed direction.

In a second test, corner radii R 1 = 0.5 mm and R 2 = 0.5 mm have been used. We can notice that the increase in the corner radius induces inclined shear bands, which

Conclusion

In this paper, a framework for practical implementation of the natural element method in the context of explicit thermo-elastoplasticity involving a strategy for dynamic adaptive refinement in this framework has been presented. The present procedure offers some outstanding advantages in large transformations including localization: (a) the C-NEM provides equivalent accuracy than the quadrilateral/hexahedral finite elements by only using the Delaunay triangulation of the current set of nodes; (b) the meshfree features of the technique allows frequent updates of the reference configuration without special attention on the relative position of the nodes; (c) the use of a stabilized conforming nodal integration simplifies the transfer of internal variables between successive updates as integration points coincide with the nodes; (d) the underlying Voronoi diagram allows refinement through unstructured scattered of nodes; (e) no size parameter is involved in the definition of the neighbors, which is essential in refinement procedure, where the nodal density is different in some parts of the domain. An example involving dynamic shear bands progression in metal has been presented to illustrate the potentiality of the proposed technique.

ΓFigure 1 .

 1 Figure 1. Construction of the Sibson shape functions.

Figure 2 .

 2 Figure 2. Support of the natural shape function related to node n i .

  iii) Perform polar decomposition of (F e ) Evaluate elastic trial logarithmic strain tensor ( e ) trial n+1 = ln[U e trial n+1 ] (v) Evaluate trial stress tensor : T trial n+1 := C e : ( e n+1 ) trial (vi) Check plastic consistency condition IF J 2 (T trial n+1 )σ n y ≤ 0 THEN : (.) n+1 = (.) trial n+1 and RETURN ELSE go to (vii) (vii) Plastic corrector (solve for Δγ)

Figure 3 .

 3 Figure 3. Representative domain around a node for the evaluation of Fi .

Figure 4 .

 4 Figure 4. Refinement procedure based on the Voronoi cells. (o): additional nodes.

Figure 5 .

 5 Figure 5. High-speed shearing: geometry of the problem.

Figure 6 .

 6 Figure 6. High-speed shearing (R 1 = 0 mm and R 2 = 0 mm): adaptive refinement during shear band formation .

Figure 7 .

 7 Figure 7. High-speed shearing (R 1 = 0 mm and R 2 = 0 mm): adaptive refinement during shear band formation and constrained Voronoi cells.

Figure 8 .Figure 9 .

 89 Figure 8. High-speed shearing (R 1 = 0 mm and R 2 = 0 mm): shear band progression (equivalent plastic strain).

Figure 10 .

 10 Figure 10. High-speed shearing (R 1 = 0.5 mm and R 2 = 0.5 mm): shear band progression (equivalent plastic strain).


Philippe Jaming
email: philippe.jaming@math.u-bordeaux.fr

Felipe Negreira
email: felipe.negreira@math.u-bordeaux.fr

A sampling theorem for functions in Besov spaces on the sphere

. The main tool is the needlet decomposition given by Narcowich et al. in [10].

I. INTRODUCTION Representation and analysis of functions (signals, images

or other data) on manifolds has become an ubiquitous task in many scientific fields. Indeed, high-dimensional data acquisition has become common but the meaningful data is often considered to belong to a high dimensional manifold in the ambient space. Among the most commonly encountered manifolds one can find the euclidean d-dimensional sphere

S d = {(x 1 , . . . , x d+1 ) ∈ R d+1 : x 2 1 + • • • + x 2 d+1 =
1}. This will be the only manifold considered in this paper. Note that the dimension here refers to the dimension of the manifold, not to that of the ambient space.

Many good representing systems (orthonormal bases, frames,...) for functions on manifolds are known (see e.g. [START_REF] Pesenson | Band-limited localized Parseval frames and Besov spaces on compact homogeneous manifolds[END_REF], [START_REF] Maggioni | Diffusion Polynomial Frames on measure metric spaces[END_REF], [START_REF] Narcowich | Decomposition of Besov and Triebel-Lizorkin spaces on the sphere[END_REF]). To properly represent a function f using such a system, one usually needs to know the set of coefficients of f in the representing system (or its dual) { f, ψ j }. However, many data acquisition systems do not provide such coefficients but rather a sample of the values of the function: {f (η j )}.

Sampling Theory precisely deals with the possibility to reconstruct a function from its samples. The most famous result in that direction is the Shannon-Whittaker-Kotelnikov Theorem which states that a band-limited function on R can be reconstructed from its regular samples. This theorem has been extended in many directions e.g. one can find a version for band-limited functions on compact homogeneous manifolds (thus including the sphere) in an article of Geller-Pesenson [START_REF] Pesenson | Band-limited localized Parseval frames and Besov spaces on compact homogeneous manifolds[END_REF]. The general idea of a sampling theorem can be stated as follows: if a function has moderate oscillation (for instance the function is band-limited, and thus has no highly oscillating Fourier component), then it is well approximated in the neighbourhood of a point by its value at that point. If those neighborhoods cover sufficiently well the space (the sampling rate is high enough), the sample allows to describe the function.

Recently, the first author and E. Malinnikova [START_REF] Jaming | An uncertainty principle and sampling inequalities in Besov spaces[END_REF] showed that, in R d a sampling theorem holds when moderate oscillation is measured in the scale of Besov spaces B s p,q . Our aim here is to show that a similar result holds on the sphere:

Main Theorem. Let 1 p ∞.
There exists finite sets of points {η 1 , . . . , η N } on S d , constants c 1 , c 2 depending only on d, and a constant M depending on d, p such that, if

f ∈ B d/p p,1 (S d ) with f B d/p p,1 M -1 f L p then c 1 f L p Nε l=1 |f (η l )| p 1/p c 2 f L p .
The actual result is more precise, see below. The proof consists in replacing the wavelet representation of Besov functions in R d by the needlet representation of Besov functions in S d as established by Narkowich-Petrushev-Ward [START_REF] Narcowich | Decomposition of Besov and Triebel-Lizorkin spaces on the sphere[END_REF] and then to adapt the computations to the spherical case.

In the next section, we introduce all material needed in the proofs so as to keep this paper self-contained. The remaining of the paper is then devoted to the proof of the main theorem. Throughout the paper, C will be a constant that changes from line to line and we add an index C d , C p , C d,p , . . . to stress the dependence on the parameters d, p, . . .. However, the constants c * , c * * , c ⋄ , κ have specific meanings stemming from the geometry of S d and are used to ease the reading of the arguments in the proof.

II. PRELIMINARIES.

In this section, we will describe a set of functions on the sphere S d that will serve as basis for various function spaces on S d . There are two choices that might be suitable for our needs. One natural choice would be the wavelet type bases that can be constructed on the sphere by following e.g. the general construction for on any space of homogeneous type (see e.g. [START_REF] Deng | Harmonic Analysis on Spaces of Homogeneous Type[END_REF]). This construction is rather evolved and only leads to a frame in L 2 (S d ). Note however that a wavelet orthonormal basis has been constructed in [START_REF] Auscher | HYTONEN Orthonormal bases of regular wavelets in spaces of homogeneous type[END_REF], but this is technically even more involved.

We therefore make the choice to use a simpler basis, namely the needlet system introduced by Narcowich, Petrushev and Ward in [START_REF] Narcowich | Decomposition of Besov and Triebel-Lizorkin spaces on the sphere[END_REF]. Let us now describe this system.

The starting point is the set of spherical harmonics. Recall that a spherical harmonic is an homogeneous polynomial that is harmonic (for the usual Laplace differential operator on R d+1 ). We denote by H ν the subspace of L 2 (S d ) of all spherical harmonics of degree ν and by P ν the orthogonal projector onto this subspace. This projection operator has an integral kernel given by a Gegenbauer polynomial, a fact that we will not need here, but allows to explicitly construct the needlet basis. We refer to [START_REF] Dai | Approximation Theory and Harmonic Analysis on Spheres and Balls[END_REF], [START_REF] Stein | Introduction to Fourier Analysis on Euclidean spaces[END_REF] for details. Anyway, this set the ground for our first definition.

Definition II.1. For each N 1 we define the kernel K N as

K N (x) := ∞ ν=0 â ν N P ν (x) (II.1)
where â is a smooth function on [0, +∞)

satisfying: i) supp â ⊂ [1/2, 2), ii) |â(t)| > c if t ∈ [3/5, 5/3], iii) â2 (t) + â2 (2t) = 1 if t ∈ [1/2, 1].
Note that condition (i) of â implies that the sum in (II.1) is finite. Moreover, ν runs over N/2 ν 2N (in particular ν 1) and therefore K N is a polynomial of degree at most 2N . While (ii) and (iii) are not necessary for (II.1) to make sense, they are needed to form the right building blocks.

These kernels have a natural action over functions defined on S d via the spherical convolution:

Definition II.2. Let Φ ∈ L ∞ [-1, 1] and f ∈ L p (S d ), 1 p ∞. We write Φ * f (ξ) := S d Φ( ξ, θ )f (θ) dσ(θ).

Next we define

Φ 0 := P 0

Φ j := K 2 j-1 for j > 0. (II.2)
We then have the following result.

Proposition II.3. For each ℓ 0 set K ℓ := ℓ j=0 Φ j . Let 1 p ∞, then for every f ∈ L p (S d ) it holds K ℓ * f -→ ℓ→+∞ f in L p (S d ). Proof. It is straightforward from the fact that ∞ ν=0 â(2 -ν t) = 1 for all t ∈ [1, +∞) and the inclusion L p (S d ) ⊂ ∞ ν=0 H ν for any 1 p ∞.
Next, in order to define our needlet system, we need to construct a particular type of partition of S d .

Lemma II.4. For any 0 < ε 1 there exists a partition R ε = {R 1 , . . . , R Nε } of S d , together with a set of points X ε = {η 1 , . . . , η Nε } ⊂ S d with the proprieties:

i) S d = R 1 ∪ • • • ∪ R Nε and the sets in R ε do not overlap ( • R l ∩ • R k = ∅ if l = k), ii) there exist a constant 0 < c * < 1 such that for every l = 1, . . . , N ε , B η l (c * ε) ⊂ R l ⊂ B η l (ε);
iii) there exist a constant c * * > 0 such that

N ε := #X ε = #R ε c * * ε -d .
The constants c * and c * * depend only on the dimension d. These sets offer a useful formula in terms of decomposition of functions.

Theorem II.6. There is a constant c ⋄ (depending only on the dimension d) such that for any N 1 and any set X ε ⊂ S d of almost uniformly ε-distributed points with ε = c ⋄ /N , there exists positive coefficients {c η } η∈Xε such that the cubature formula

S d p(ξ)dξ ∼ η∈Xε c η p(η)
is exact for all spherical polynomials p(ξ) of degree N .

In addition for all η, c 1 N -d c η c 2 N -d with the equivalence constants c 1 , c 2 only depending on d.

For the construction of the needlet system, we are going to denote for each j 0, ε j := c ⋄ /2 j+2 , X j a fixed set of ε j -distributed points and {c η } η∈Xj its associated family of positive coefficients. Also, we set X := ∞ j=0 X j .

Then, from (II.2), we define our needlet system for each j 0 as

ϕ η (ξ) := √ c η Φ j (η • ξ), η ∈ X j . (II.3)
Adding all up we have the following decomposition result.

Theorem II.7.

Let f ∈ L p (S d ), 1 p ∞. Then f = η∈X f, ϕ η ϕ η in L p (S d ).
(II.4)

Sketch of proof.

Using the identity f = ∞ j=0 Φ j * Φ j * f , we decompose Φ j * Φ j (ξ • •) for each j with Theorem II.6 and then we sum all over X = ∞ j=0 X j to get (II.4). The details of the proof can be found in [START_REF] Narcowich | Decomposition of Besov and Triebel-Lizorkin spaces on the sphere[END_REF]Proposition 3.2].

The main difference that needlets have with respect to wavelets is that their localization is only measured in terms of a decay property and not by their support. More precisely, they satisfy the following size condition: for every k ≥ 0, there exist a constant c k such that for every j 0 and η ∈ X j it holds

|ϕ η (ξ)| c k 2 jd/2 (1 + 2 j d(ξ, η)) k , ∀ξ ∈ S d .
(II.5)

And they also have a smoothness condition: there exist another constant κ such that for every j ≥ 0, if d(ξ, θ) κ2 -j and η ∈ X j then

|ϕ η (ξ) -ϕ η (θ)| c k 2 jd/2 2 j d(ξ, θ) (1 + 2 j d(ξ, η)) k . (II.6)
We can now define the Besov spaces in S d . Again, we follow [START_REF] Narcowich | Decomposition of Besov and Triebel-Lizorkin spaces on the sphere[END_REF].

Definition II.8. Given 0 < p, q ∞ and s ∈ R, the Besov space B s p,q (S d ) is defined as the set of all f ∈ S ′ such that the norm

f B s p,q :=   ∞ j=0 2 js Φ j * f L p q   1/q
is finite. As usual, the L p , ℓ q norms are replaced by the supnorms when p = ∞ or q = ∞.

From this definition it is easy to see that B s ′ p,q ⊂ B s p,q when s ′ s, and B s p,q ′ ⊂ B s p,q when q ′ q. However, for our purposes, it will be more useful to characterize this space through the needlet decomposition

f B s p,q :=    ∞ j=0   2 j(s+d/2-d/p)   η∈Xj | f, ϕ η | p   1/p    q    1/q
It is important to remark that these two definitions indeed coincide, for a proof of this equivalence see [

(M ε α ) -1 f L p then 1 2ε d/p f L p Nε l=1 |f (ζ l )| p 1/p 3C 2ε d/p f L p .
Proof. From (II.4) we may write f ∈ B s p,1 (S d ) as

f = ∞ j=0 η∈Xj f, ϕ η ϕ η
where ϕ η are the needlets defined in (II.3). Let {R 1 , . . . , R Nε } be the partition of S d associated to {ζ 1 , . . . , ζ Nε }, take an arbitrary l ∈ {1, . . . , N ε } and consider ξ ∈ R l . We then have

|f (ξ) -f (ζ l )| ∞ j=0 η∈Xj | f, ϕ η ||ϕ η (ξ) -ϕ η (ζ l )| ∞ j=0   η∈Xj |ϕ η (ξ) -ϕ η (ζ l )| p ′   1/p ′ E j (f ) (III.1)
where for each j 0

E j (f ) :=   η∈Xj | f, ϕ η | p   1/p .
In order to get the right bound for (III.1), we are going to use two estimates which we will prove after showing that they imply the result.

Claim 1. For every j 0 and every ξ ∈ S d ,

  η∈Xj |ϕ η (ξ)| p ′   1/p ′ C p 2 jd/2
where C p depends only on d and p.

In particular this implies that

  η∈Xj |ϕ η (ξ) -ϕ η (ζ l )| p ′   1/p ′ 2C p 2 jd/2 C p 2 -j0(s-d/p) 2 j(s+d/2-d/p)
where 2 j ≥ 2 j0 with j 0 to be chosen later.

Claim 2. For every j such that 2 j κε -1 and ξ

∈ B ζ l (ε),   η∈Xj |ϕ η (ξ) -ϕ η (ζ l )| p ′   1/p ′ C p 2 jd/2 2 j d(ξ, ζ l )
where C p depends only on d and p.

Note that since R l ⊂ B ζ l (ε) by propriety (ii) in Lemma II.4, the condition ξ ∈ B ζ l (ε) is always met when ξ ∈ R l .
Assuming the claims and chosing j 0 so that 2 j0 κε α ε -1 κε -1 , (III.1) reads

|f (ξ) -f (ζ l )| C p j0 j=0 2 jd/2 2 j d(ξ, ζ l )E j (f ) + C p 2 -j0(s-d/p) ∞ j=j0+1 2 j(s+d/2-d/p) E j (f ).
Next, taking the L p -norm over R l and using the triangular inequality, we get

R l |f (ξ) -f (ζ l )| p dσ(ξ) 1/p C p j0 j=0 2 jd/2 2 j R l d(ξ, ζ l ) p dσ(ξ) 1/p E j (f ) + C p 2 -j0(s-d/p) σ(R l ) 1/p ∞ j=j0+1 2 j(s+d/2-d/p) E j (f ) C p j0 j=0 2 jd/2 2 j c d ε d/p+1 E j (f ) + C p 2 -j0(s-d/p) ε d/p ∞ j=j0+1 2 j(s+d/2-d/p) E j (f )
where we use again the fact that R l ⊂ B η l (ε). We now take 2 j0 ≈ ε α-1 so that 2 -j0(s-d/p) Cε α and thus

R l |f (ξ) -f (ζ l )| p dσ(ξ) 1/p C p 2 j0 εε d/p j0 j=0 2 j(s+d/2-d/p) E j (f ) + C p 2 -j0(s-d/p) ε d/p ∞ j=j0+1 2 j(s+d/2-d/p) E j (f ) C p ε α ε d/p f B s p,1
which combined with the triangular inequality gives us

R l |f (ξ)| p dσ(ξ) 1/p -|f (ζ l )|σ(R l ) 1/p C p ε α ε d/p f B s p,1 .
So, by taking ℓ p -norms and using the triangular inequality one last time, we conclude that

f L p - Nε l=1 |f (ζ l )| p σ(R l ) 1/p = R l |f (ξ)| p dσ(ξ) 1/p ℓ p Nε -|f (ζ l )|σ(R l ) 1/p ℓ p Nε R l |f (ξ)| p dσ(ξ) 1/p -|f (ζ l )|σ(R l ) 1/p ℓ p Nε C p ε α ε d/p f B s p,1 N 1/p ε C p (c * * ) 1/p ε α f B s p,1
where c * * is the constant in property (iii) of Lemma II.4.

Then we have

Nε l=1 |f (ζ l )| p σ(R l ) 1/p C p (c * * ) 1/p ε α f B s p,1 + f L p and Nε l=1 |f (ζ l )| p σ(R l ) 1/p f L p -C p (c * * ) 1/p ε α f B s p,1 .
Therefore, choosing M = 2C p (c * * ) 1/p , we have that

1 2 f L p Nε l=1 |f (ζ l )| p σ(R l ) 1/p 3 2 f L p whenever f B s p,1 (M ε α ) -1 f L p .
Which proves the theorem, since by the propriety (ii) in Lemma II.4 we know that

B η (c * ε) ⊂ R l ⊂ B η (ε) and thus c * ε d σ(R l ) ε d .
Proof of the Claims. Let us first prove Claim 1. Using the size condition (II.5) for ϕ η with k = d + 1 we have that

  η∈Xj |ϕ η (ξ)| p ′   1/p ′ c d+1 2 jd/2   η∈Xj 1 (1 + 2 j d(ξ, η)) (d+1)p ′   1/p ′ . (III.2)
So, in order to get the bound in Claim 1, it is enough to prove that the sum on the right hand side is bounded by a constant depending only on d and p. To that end, define for each η ∈ X j , G η := B η (c * 2 -j ) the ball in (ii) of Lemma II.4. The triangular inequality shows that

1 + 2 j d(ξ, η) (1 -c * ) max θ∈Gη (1 + 2 j d(ξ, θ)).
And then, as {G η } η∈Xj is a family of disjoint sets of measure (c * ) d 2 -jd , summing over all X j gives us

η∈Xj 1 (1 + 2 j d(ξ, η)) (d+1)p ′ η∈Xj (1 -c * ) -(d+1)p ′ σ(G η ) Gη 1 (1 + 2 j d(ξ, θ)) (d+1)p ′ dσ(θ) (1 -c * ) -(d+1)p ′ S d 2 jd (1 + 2 j d(ξ, θ)) (d+1)p ′ dσ(θ) c p (1 -c * ) -(d+1)p ′
by standard calculus. This completes the proof of Claim 1. The proof of Claim 2 is exactly the same using (II.6) instead of (II.5) in (III.2).

As the needlets ϕ η are continuous functions (moreover they are polynomials) the same computations above and the fact the S d is compact show that every f ∈ B s p,1 (S d ) is continuous. Also, another corollary that can be drawn from the previous calculations is the following result: holds true for all f ∈ B s p,1 (S d ).

IV. CONCLUSION AND FUTURE WORK

In this paper, we have used the needlet decomposition of Besov spaces to establish a sampling theorem in B s p,1 when the smootheness index is large enough. This opens the path to future work in two directions:

-establish efficient algorithms for reconstructions of function on a sphere from their samples;

-extend this work to more general compact Riemannian manifolds for which needlet like decompositions have been established in [START_REF] Pesenson | Band-limited localized Parseval frames and Besov spaces on compact homogeneous manifolds[END_REF], [START_REF] Maggioni | Diffusion Polynomial Frames on measure metric spaces[END_REF] and then to more general homogeneous spaces using wavelet decompositions like those given in [START_REF] Deng | Harmonic Analysis on Spaces of Homogeneous Type[END_REF]. Note that some results in this direction have already been obtained by various authors. For instance H. Führ and K. Gröchenig [START_REF] Öchenig | Sampling theorems on locally compact groups from oscillation estimates[END_REF] obtained sampling theorems under oscillation conditions on locally compact groups while I.Z. Pesenson [START_REF] Pesenson | Paley-Wiener approximations and multiscale approximations in Sobolev and Besov spaces on manifolds[END_REF] established some links between sampling and smallness of Besov norms on compact Riemannian manifolds.

  Throughout this paper, we will adopt the following notation: on R d+1 , |•| and •, • are the usual Euclidean norm and scalar product while S d = {x ∈ R d+1 : |x| = 1} is the unit sphere. We endow S d with its usual geodesic distance d(ξ, θ) = 1 π arccos ξ, θ and with σ the Lebesgue measure. The balls on S d are denote by B η (r) := {ξ ∈ S d : d(ξ, η) r} where η ∈ S d and r ≤ 1.

For a proof of this lemma see [ 9 ,

 9 Proposition 3.2]. It is a particular case of a more general construction of Hytönen and Kairema [6, Theorem 2.1]. Definition II.5. A set X ε ⊂ S d which along with the associated partition R ε of S d has the proprieties given in Lemma II.4 will be be called a set of almost uniformly εdistributed points on S d .

Corollary III. 2 .

 2 Let 1 p ∞, s ≥ d/p and α = s-d/p 1+s-d/p . Fix ε > 0 and let {ζ 1 , . . . , ζ Nε } be a set of almost uniformly ε-distributed points and let {R 1 , . . . , R Nε } be the associated partition of S d . Then there exists a constant C p depending only on d and p such thatS d f (ξ) -Nε l=1 f (ζ l )1 R l (ξ) p dσ(ξ) C p ε α f B s p,1


Philippe Jaming
email: philippe.jaming@math.u-bordeaux.fr

Phase retrieval for solutions of the Shrödinger and Helmholtz equations

In this paper we give an overview of our recent results on the phase retrieval problem for solutions of two partial differential equations: the free Shrödinger equation and the Helmholtz equation. In particular, we show that structured signals often imply uniqueness in the phase retrieval problem.

Phase retrieval;Shrödinger equation; Helmholtz equation

I. INTRODUCTION

The phase retrieval problem consists in reconstructing a function from its modulus or the modulus of some transform of it (frame coefficient, Fourier transform,...) and some structural information on the function (e.g. to be compactly supported). Such a problem occurs in many scientific fields: microscopy, holography, crystallography, neutron radiography, optical coherence tomography, optical design, radar signal processing, quantum mechanics to name a few. We refer to the books [Hu], [St], the review articles [KST], [Mi], [Fi], [LBL] and to the introduction of our previous paper [Ja] for descriptions of various instances of this problems, for some solutions to it (both theoretical and numerical) and for further references.

The problem can be split into two main questions:

-design algorithms for the reconstruction of at least one solution.

-obtain uniqueness results.

After having long been ignored by mathematicians, recent progress on the algorithmic aspect of the problem [CSV], [WdAM] has triggered a lot of attention to this problem. While the design of numerical algorithms that allow the reconstruction of one solution is of course essential, the task can only be complete once one is certain to reconstruct all solutions of interest. This is generally not possible as long as uniqueness is not guarantied and plainly justifies the second part of the problem. Uniqueness is also usefull in order to have stability results in presence of additive noise. In this paper, we will only deal with the uniqueness aspect of the problem. More precisely, the phase retrieval problem is extremely common in optical sciences. Among other reasons, this is due to the lack of sensitivity of optical measurement instruments to phase. It turns out that optical signals are solution of Partial Differential Equations and our aim is to show that this information can be of some use in the phase retrieval problem. We will concentrate on two equations.

The first equation we consider is the free Shrödinger equation on R d with initial condition

u 0 ∈ L 2 (R d ) :    i∂ t u + 1 4π ∆ 2 x u = 0 u(x, 0) = u 0 (x) (I.1)
where ∆ x = j ∂ 2 xj is the usual Laplace operator. In this case, we ask the following: An entire (defined on the whole of R d ) solution u of (I.3) can be expressed in polar coordinates as a series

Problem 1 (Phase Retrieval Problem for the Free Shrödinger Equation). Let u 0 , v 0 ∈ L 2 (R d
u(rθ) ∼ (2π) 1/2 r (d-2)/2 ∞ m=0 N (m) j=1 a m,j J ν(m) (r)Y j m (θ) (I.4)
where J ν are the Bessel functions, ν(m) = m + (d -2)/2 and {Y j m } j=1,...,N (m) is a basis for the spherical harmonics of degree m in R d .

The far-field-pattern f is defined to be the spherical part of r (d-1)/2 u(rθ) as r → ∞, i.e., for θ ∈ S d-1 (the unit sphere of R d ),

f (θ) ∼ ∞ m=0 N (m) j=1 a m,j (f )Y j m (θ). (I.5)
A solution u of (I.4) is said to be a Herglotz wave function if its far-field-pattern is a function f ∈ L 2 (S d-1 ) and we then write u = u f . We will denote the space of Herglotz wave functions by W 2 (R d ). The problem we then address is the following:

Problem 2 (Phase Retrieval Problem for the Helmholtz Equation).

Let u f , u g ∈ W 2 (R d ) be such that |u g | = |u f |. Does this imply that either u g = cu f or u g = cū f for some c ∈ C with |c| = 1. What if f is in some subspace of L 2 (S d-1
), or if both f, g belong to some subspace L 2 (S d-1 ).

In this paper, we will report on our work [Ja], [AJ], [JPE] to show that, under some structural restriction on the solution, uniqueness may be achieved in the phase retrieval problem.

II. THE FREE SHR ÖDINGER EQUATION

A. Reformulation of the problem

Let us recall that the solution of the Free Shrödinger Equation (I.1) with initial data

u 0 ∈ L 2 (R) has solution u(x, t) = R e -iπ|ξ| 2 t+2iπ x,ξ u 0 (ξ) dξ = F -1 e -iπξ 2 t u 0 (x)
where F -1 is the inverse Fourier transform. Let us introduce the Fractional Fourier Transform (FrFT

): let α ∈ R \ πZ, let c α = exp i 2 α -π 2 | sin α| . For u ∈ L 1 (R d ) and α / ∈ πZ, define F α u(ξ) = c α e -iπ|ξ| 2 cot α × × R u(t)e -iπt 2 cot α e -2iπtξ/ sin α dt.
A straightforward computation then shows that, for α ∈ (-π/2, π/2),

F α u 0 (ξ) = ie iα/2 | cos α| d e -iπ|ξ| 2 cot α u(ξ/ cos α, cot α).
Next, let us now recall that the ambiguity function of u ∈ L 2 (R) is defined by

A(u)(x, y) = R u t + x 2 u t - x 2 e -2iπty dt.
The link between the FrFT and the ambiguity function comes from the following property:

A(F α u)(x, y) = A(u)(x cos α -y sin α, x sin α + y cos α).
In particular,

A(u)(-y sin α, y cos α) = A(F α u)(0, y) = F [|F α u| 2 ](y). (II.1)
It follows that the knowledge of |u(x, t)| is equivalent to the knowledge of A(u) on the line R(1, t). In other words, the Phase Retrieval Problem 1 has been changed into a sampling problem for the ambiguity function. Exploiting this link, we have been able to prove the following:

Theorem II.1 (Jaming [Ja]). Let u 0 , v 0 ∈ L 2 (R) and u, v be the corresponding solutions of (I.1). Let τ ⊂ R be a set of times and assume that

|u(x, t)| = |v(x, t)| for every x ∈ R d and t ∈ τ. 1) If τ = R (or is dense in R). Then there exists c ∈ C with |c| = 1 such that v = cu.
2) If τ ⊂ R has an accumulation point (in particular if τ has positive measure) and if u 0 , v 0 are both compactly supported, then there exists c

∈ C with |c| = 1 such that v = cu. 3) Let a > 0 if u, v ∈ L 2 (R) have compact support included in [-a, a] and if τ = k a 2 , k ∈ Z then there exists c ∈ C with |c| = 1 such that v = cu. Sketch of proof. If |u(x, t)| = |v(x, t)| for every x, t then A(u) = A(v)
and it is well known that this implies u = v.

The two other items are based on the observation that

A(u)(x, •) is the Fourier transform of the function ϕ x (t) := u t + x 2 u t -x 2 . But, if u is compactly supported, so is ϕ x , thus A(u) is an entire function in the second variable thus if A(u)(x, •), A(v)(x,
•) coincide on a large set, they coincide everywhere.

Further, if u and v are supported in intervals, so is ϕ x and we may simply apply Shannon's sampling theorem. This theorem still requires the knowledge of the solution at infinitely many times. It turns out that if we ask u, v to be "structured", then less is needed:

Theorem II.2 (Jaming [Ja]). Let u 0 , v 0 ∈ L 2 (R) and u, v be the corresponding solutions of (I.1). Let τ ⊂ R be a set of times and assume that

|u(x, t)| = |v(x, t)|
for every x ∈ R d and t ∈ τ.

1) Assume u 0 , v 0 are Hermite functions (of the form P γ where P is a polynomial and γ

(x) = e -πx 2 ), τ = {cot α} with α ∈ R \ Qπ. Then v = cu with |c| = 1.
2) Assume u 0 , v 0 are rectangular pulse trains (of the form

f inite c j 1 [j,j+η) with 0 < η ≤ 1
2 and c j a finite sequence of complex numbers) and τ = {t 0 }, t 0 > 0, then v = cu with |c| = 1.

3) Assume u 0 , v 0 are of the form f inite c i e -π(t-ai) 2 , c i ∈ C and a i ∈ R, and τ = {t 0 }, t 0 > 0, then v = cu with |c| = 1.

Sketch of proof.

The proofs of these three facts are algebraic in nature and require to expand u 0 , v 0 in some representing system and to identify coefficients.

For Hermite functions, the problem is best solved for the FrFT and we use the Hermite basis and write u(t)

= N j=0 c j H j , v(t) = M j=0
c j H j . Then, as is well known,

F θ [u](t) = N k=0 c k e -ikθ H k (t)
and a similar expression holds for

F θ [v]. Then |v| 2 = |u| 2 and |F θ [v]| 2 = |F θ [u]| 2 is equivalent to M j,k=1 d j d k H j (t)H k (t) = N j,k=1 c j c k H j (t)H k (t)
and M j,k=1

d j d k e i(k-j)β H j (t)H k (t) = N j,k=1 c j c k e i(k-j)β H j (t)H k (t)
An induction starting from the highest order term then leads to the result. For the pulse trains, the ambiguity function reformulation is the most convenient. We first introduce the discrete ambiguity function of a sequence a:

A(a)(k, y) = j∈Z a j a j-k e 2iπjy .
Then we write u(t) = a j 1 [0,η) (tj) where a j is a finite sequence. A simple commutation shows that, for x ∈

j - 1 2 , j + 1 2 , A(u)(x, y) = e iπjy A(a)(j, y)A(1 [0,η) )(x -j, y)
and a similar expression holds for A(v). The part A(1 [0,η) )(xj, y) being common to A(u) and A(v), the problem reduces to a problem on the discrete part which is again solve by an induction procedure.

The last case is more elaborate and we refer to [Ja].

This theorem mislead us to conjecture that three different times might be sufficient to obtain uniqueness. This was disproved in [CHST] who noticed that as soon as there exists a triplet for which non-uniqueness could be achieved, then no triplet would work and then elaborated on a counter-example given in [Ja] that showed that the condition on the time in the case of Hermite function is necessary. Let us here give a different proof that appeals to an argument of A.E.J.M Janssen [Jan] and an oblique marginal formula for the Zak transform

Z[f ](x, ξ) = k∈Z f (x + k)e -2iπkξ .
which is of independent interest.

Theorem II.3 (Andreys-Jaming [AJ]). Let p ∈ Z \ {0}, q ∈ N be relatively prime integers and let α be defined by cot α = p q . For n ∈ Z, and x ∈ [0, 1], let c n,p,q = 1 q q-1 k=0 (-1) kp e iπ p q k 2 e -2iπ nk q and, for x ∈ [0, 1], n ∈ Z, write ξ n,p,q (x) = p q x + 1 2 p + n q and A p,q (x) = {n ∈ Z| 0 ≤ ξ n,p,q (x) ≤ 1}. Then, for every f ∈ L 2 (R),

F α f (ξ) = c α e -iπ p q ξ 2 1 0 e -2iπ ξ sin α x-iπ p q x 2 × n∈Ap,q(x) c n,p,q Zf x, ξ sin α + ξ n,p,q (x) dx. (II.2)
The identity has to be taken in the L 2 (R) sense.

This theorem is a (rigourous) elaboration on the formal computation that allows to compute the Zak transform of a chirp. From this, it is easy to obtain the following:

Corollary II.4 (Andreys-Jaming [AJ]). Let r 1 , . . . , r m ∈ Q be m different rational numbers and N ∈ N. For j = 1, . . . , m let α j ∈ [0, π) be defined via cot α j = r j . Then there exists

f 1 , . . . , f N ∈ L 2 (R) such that, for j = 1, . . . , m and k = ℓ = 1, . . . , N , supp F αj f k ∩ supp F αj f ℓ = ∅. (II.3)
In particular, for every j = 1, . . . , m and for every

c 1 , . . . , c N ∈ C with |c k | = 1, F αj N k=1 c k f k = F αj N k=1 f k .
(II.4)

Taking u 0 = N k=1 f k , v 0 = N k=1 c k f k
and solving the Shrödinger equation with these initial conditions, we obtain functions u(t, x), v(t, x) that have modulus coinciding at m different times but still being different.

Corollary II.4 is established as follows: if one takes a function F ∈ L 2 supported in a (union) small disc in [0, 1] × [0, 1], then this function is the Zak transform of some function f ∈ L 2 (R), F = Zf and the support of F α [f ] is then obtained from (II.2) by looking at all oblique lines (x, ξ sin α + ξ n,p,q (x)) appearing in that formula that intersect the disc (see Figure 1).

III. THE HELMHOLTZ EQUATION

The aim of this section is to show that solutions of the Helmholtz equation are often uniquely determined by their modulus. The main result is the following:

Theorem III.1 (Jaming & Pérez-Esteva [JPE]). Let u f , u g ∈ W 2 (R d ) be such that |u g | = |u f |.
Assume one of the following holds:

f and g are real valued; -f has non-zero mean; -the dimension is d = 2; ¨¨¨¨¨¨¨¨¨¨¨z z z ~u J J Fig. 1. In this figure p = 2, q = 1 and the slope is

p q = 1 2 . The support of Z[f ]
is the union of the black discs. Then the set J is such that, for each η ∈ J, the (periodic) line starting at (0, η) with slope q/p = 1/2 (tan α = q/p) intersects at least once the support of Zf (one such line is drawn dotted). Here J = [0, 1/4] ∪ [1/2, 3/4] Finally, the support of

Fα[f ] is then included in the set of all ξ such that ξ/ sin α + p/2 ∈ J + Z = k∈Z [k/2, (k + 1/2)/2]. As p = 2, q = 1, sin α = 2/ √ 5, thus supp Fα[f ] ⊂ k∈Z [k/ √ 5, (k + 1/2)/ √ 5].
d ≥ 3 and f, g are zonal functions. Then either

u g = cu f or u b = cu f for some c ∈ C with |c| = 1.
Recall that a zonal function is a function of the form ϕ( x, x 0 ) for some x 0 ∈ S d-1 . Sketch of proof. The case f, g being real is of a slightly different nature to the other statements. In this case, u f , u g are also real and

|u f | = |u g | is equivalent to u 2 f = u 2 g thus (u f -u g )(u f +u g ) = 0.
It follows that at least one of u f = u g or u f = -u g holds on a set of positive measure. But they are also analytic, so that if such a relation holds on a set of positive measure, it holds everywhere.

Next notice that if

u f ∈ W 2 (R d ), then (I.4) implies |u f (rθ)| 2 = 2π r d-2 ∞ m,n=0 c m,n (f )J m+ d-2 2 (r)J n+ d-2 2 (r) where c m,n (f ) = N (m) j=1 N (n) k=1 a m,j (f )a n,k (f )Y j m (θ)Y k n (θ).
The key observation is that, thanks to the properties of the Bessel functions, |u f | = |u g | can be reformulated in terms of the c m,n 's only:

Lemma III.2. Let u f , u g ∈ W 2 (R d ). Then |u f | = |u g | if and only if for every 0 ≤ m ≤ n ℜ c m,n-m (f ) = ℜ c m,n-m (g) .
(III.1)

Then a 0,1 (f ) is the mean of f and |a 0,1 (f )| 2 = c 0,0 (f ) = c 0,0 (g) = |a 0,1 (g)| 2 . Up to replacing f, g by multiples of them, we can assume that a 0,1 (f ) = a 0,1 (g) is real (and non zero). Then ℜ c 0,n (f ) = ℜ c 0,n (g) for every n.

But, if we chose the basis of spherical harmonics to be real valued, this implies that ℜ a j,k (f ) = ℜ a j,k (f ) for every j, k thus ℜ(f ) = ℜ(g) and ℜ(u f ) = ℜ(u g ). As One then uses the other relations to show that κ m does not depend on m and that the choice between both alternatives also does not depend on m.

  ) and let u and v be the solutions of the Free Shrödinger equation (I.1) with initial value u 0 and v 0 . Let τ = {t i } i∈I ⊂ [0, +∞) be a set of times (finite or not) at which one measures |u| and |v|. Assume that |u(x, t i )| = |v(x, t i )| for all i ∈ I. (I.2) i) Does this imply that u = cv for some constant c ∈ C, |c| = 1? ii) If we restrict u ∈ D for some set D ⊂ L 2 (R d ) does (I.2) imply that u = cv for some constant c ∈ C, |c| = 1? iii) If we further restrict both u, v ∈ D for some set D ⊂ L 2 (R d ) does (I.2) imply that u = cv for some constant c ∈ C, |c| = 1? The second equation we consider is the Helmholtz equation on R d ∆u + u = 0. (I.3)

  |u g | 2 = |u f | 2 we deduce that |Im (u g )| 2 = |Im (u f )| 2 . But Im u g = u Im (g) , Im u f = u Im (f ) are real Herglotz functions. The real valued case then implies that either Im u g = Im u f or Im u g = -Im u f , that is g = f or g = f .The case of zonal functions is of similar flavor when d ≥ 3. The d = 2 case uses the fact that ℜc m,m (f ) = ℜc m,m (g) reduces to | f (m)e imθ + f (-m)e -imθ | = | g(m)e imθ + g(-m)e -imθ | so that, for each m, there exists κ m ∈ C with |κ m | = 1 such that -either f (m)e imθ + f (-m)e -imθ = κ m g(m)e imθ + g(-m)e -imθ or f (m)e imθ + f (-m)e -imθ = κ m g(m)e imθ + g(-m)e -imθ .


come

Evaluation and correction of angular anisotropic effects in multidate SPOT and Thematic Mapper data

Open Archive TOULOUSE Archive Ouverte (OATAO)

OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible. SPOT and Thematic Mapper data. (1993) Remote Sensing of Environment, Any correspondence concerning this service should be sent to the repository administrator: staff-oatao@listes-diff.inp-toulouse.fr

  This is an author-deposited version published in: http://oatao.univ-toulouse.fr/ Eprints ID : 4430

	To link to this article : DOI:10.1016/0034-4257(93)90112-B
	URL : http://dx.doi.org/10.1016/0034-4257(93)90112-B
	To cite this version : Muller, Etienne Evaluation and correction of
	angular anisotropic effects in multidate


Elli Zavou

Ordonnancement en ligne pour les machines parall èles

Keywords: ordonnancement des tâches, algorithmes en ligne, tolérance aux pannes, augmentation des ressources

L'exécution stable des tâches dans les machines parallèles est très importante. Par leur nature dynamique, ces systèmes doivent faire face à un défi de taille : ils doivent pouvoir répondre aux continuelles requêtes des utilisateurs, qui peuvent requérir des traitements différenciés. De plus, ces requêtes peuvent subir des erreurs imprévisibles, produites soit par un équipement malintentionnée, soit par un taux d'arrivées trop élevé. La consommation d'énergie induite peut également s'avérer importante, ce qui présente un autre défi. Dans ce papier nous considérons ces deux défis et conduisons une analyse compétitive au pire cas des performances d'algorithmes déterministes en ligne. Nous supposons également une sorte d'augmentation des ressources, d'accélération de la machine, qui caractérise la consommation d'énergie du système. Pour les mesures de performances, nous utilisons la charge complète, la charge en attente, ainsi que le ratio de la latence par rapport aux tâches réalisées. Nous montrons qu'il existe un seuil d'accélération en dessous duquel aucune compétitivité ne peut être atteinte par les algorithmes déterministes, même dans le cas d'une seule machine, et au delà duquel nous analysons les performances des algorithmes les plus utilisés et proposons de nouveaux algorithmes que nous démontrerons comme étant optimaux.

Introduction

Motivation. The fast development of computing systems as well as the increase of computationallyintensive demands from the users, are just two of the triggers that have led to the development of multicorebased parallel machines, internet-based computing platforms and co-operational distributed systems. Apart from the dynamic requests from the users (tasks arrivals) that may have different computational requirements (task sizes), these systems often suffer from unpredictable processor failures (machine crashes and restarts), either malicious or due to overload. In this work we choose to consider speed augmentation in order to overcome these challenges, increasing the computational power of the machines to a speedup s ≥ 1. Under this speedup, the machines execute a task s times faster than its baseline execution time (size). It affects however, the power consumption of the system. We are hence interested in the trade-off between the speedup and the guaranteed performance of the scheduling algorithms. Our main goal is to achieve reliable and stable computations while keeping the energy consumption of the system to a minimum.

Model. We consider a system of m homogeneous, fault-prone machines, with unique ids from the set [m] = {1, 2, . . . , m}. We assume that they have access to a shared object, called Repository, which represents the interface of the system with the users. The users submit tasks to the repository, the machines decide in a parallel way which tasks they schedule and after each successful execution they notify the repository of their completion. The users are then notified accordingly. A task τ has an arrival time α(τ) and a size c(τ) ∈ [c min , c max ], where c min and c max the smallest and largest size possible, and ρ = c max /c min the size ratio. Note that the size of a task represents the time it needs to be executed by a machine running with no speedup, i.e., s = 1, and that task executions are atomic, meaning that preemption and migration are not allowed. These task arrivals are represented by an infinite arrival pattern A, and the number of task sizes considered is denoted by parameter k. Furthermore, the machines may suffer from crashes and restarts, which are represented by an infinite failure pattern E. The task that is being executed by the machine that is crashed at a certain instant is not completed and has to be re-scheduled eventually. An adversarial entity is assumed to have control of both A and E in such a way that it gives worst-case scenarios, however assuming that at every point in time at least one machine will be available (not crashed). Due to the three main parameters of the model : number of machines, m, amount of speedup, s, and number of different task sizes, k, we denote it by M m, s, k .

We focus on three evaluation metrics embodying the machine utilization with completed load, the buffering with pending load, and the user fairness of the system with latency. The completed load is defined as C s t (ALG, A, E) = ∑ w∈N s t (ALG,A,E) c(w), the pending load as P s t (ALG, A, E) = ∑ w∈Q s t (ALG,A,E) c(w) and the latency as

L s t (ALG, A, E) = max{ f (w) -α(w) : ∀w ∈ N s t (ALG, A, E),t -α(w) : ∀w ∈ Q s t (ALG, A, E)},
where N and Q are the sets of completed and pending tasks respectively and f (w) is the completion time of task w. Note that, computing the optimal schedule for the three measures offline (knowing a priori patterns A and E) is an NP-hard problem [START_REF] Zavou | Online Scheduling in Fault-prone Systems : Performance Optimization and Energy Efficiency[END_REF].

Since the scheduling decisions are to be made in a continuous manner and without knowledge of any future information, we see the problem as an online scheduling problem [START_REF] Pruhs | Online scheduling. handbook of scheduling : Algorithms, models, and performance analysis[END_REF] and perform asymptotic competitive analysis [START_REF] Borodin | Online Computation and Competitive Analysis[END_REF][START_REF] Van Stee | On-line scheduling and bin packing[END_REF] to evaluate the performance of deterministic and work-conserving online algorithms, under worst-case scenarios. The asymptotic performance ratio that corresponds to the analysis, is the long-term competitive ratio over the sets of arrival and error patterns, A and E, against any algorithm in the set of algorithms X that solves the scheduling problem. Note here that we only consider the combinations of patterns for which the completed load of any algorithm X ∈ X goes to infinity, i.e., lim t→∞ C s t (X, A, E) = ∞. Hence, the three performance measures are defined as :

Completed Load : C s (ALG, A,E) = inf A∈A,E∈E,X∈X lim t→∞ C s t (ALG,A,E) C 1 t (X,A,E) .
Pending Load :

P s (ALG, A,E) = sup A∈A,E∈E,X∈X lim t→∞ P s t (ALG,A,E) P 1 t (X,A,E)
.

Latency : L s (ALG, A,E) = sup A∈A,E∈E,X∈X lim t→∞ L s t (ALG,A,E) L 1 t (X,A,E) .
Contributions. This work presents some of the most important results of my thesis [START_REF] Zavou | Online Scheduling in Fault-prone Systems : Performance Optimization and Energy Efficiency[END_REF], focusing on deterministic, work-conserving algorithms. It shows some general results that hold for all online deterministic scheduling algorithms and then focuses on the performance of four popular scheduling algorithms, analyzing their fault-tolerant properties under the worst-case scenarios created by the adversarial entity. The algorithms are : the Longest In System (LIS), which schedules the task that has been waiting the longest in the repository, the Shortest In System (SIS), which schedules the task that has arrived latest, the Largest Processing Time (LPT), which schedules the pending task of the biggest size, and the Shortest Processing Time (SPT), which schedules the pending task of the smallest size. Finally, some alternative algorithms are proposed, that achieve optimal competitiveness in the specified models.

Results

Table 1 summarizes the results obtained for any deterministic and work-conserving scheduling algorithm, as well as some more specific results for the four widely-used algorithms, in the case of a single machine.It gives an insight to the reader for the challenges of online scheduling algorithms even in the simplest model, thus providing a first idea of the limitations they will have in the case of m parallel machines. Each row represents the performance of an algorithm in the specified model with respect to the three evaluation metrics, and provides the reference of the published paper where their detailed analysis can be found.

The first group of results is about all deterministic and work-conserving algorithms, where one can clearly see the limitations even in the case of one machine. In the first row, by assuming no speedup and an infinite amount of task sizes available, we show that no such algorithm can achieve any competitiveness. In the second row, we bound the amount of task sizes available to two, and we are able to show that algorithms may only achieve up to 1/2-completed-load competitiveness. In the third row, we allow some speedup instead of bounding the number of task sizes, and show that no algorithm can be 1-completed-load competitive (cannot be optimal) if the speedup is less than the minimum of ρ and 1 + γ/ρ, where γ = max{ ρ-s s-1 , 0}.

Ordonnancement en ligne pour les machines parallèles

Alg.

Model Completed

Pending Latency, L Reference Load, C Load, P Determ.

M 1, 1, ∞ 0 ∞ ∞ [FAGKZ16, KWZ15] M 1, 1, 2 ≤ ρ ρ+ ρ ≈ 1 2 ∞ ∞ [FAGK + 15, FAGKZ15] M 1, s < min{ρ, 1 + γ/ρ}, ∞ < 1 ∞ ∞ [FAGKZ16, KWZ15] LIS M 1, s < ρ, 2 0 ∞ ∞ [FAGKZ16] M 1, s ∈ [ρ, 1 + 1/ρ), ∞ [1/ρ, 1 2 + 1 2ρ ] [ 1+ρ 2 , ρ] (0, 1] [FAGKZ16] M 1, s ∈ [max{ρ, 1 + 1 ρ }, 2), ∞ [1/ρ, s/2] [ s 2(s-1) , ρ] (0, 1] [FAGKZ16] M 1, s ≥ max{ρ, 2}, ∞ 1 1 (0, 1] [FAGKZ16] SIS M 1, s < ρ, 2 0 ∞ ∞ [FAGKZ16] M 1, s ∈ [ρ, 1 + 1/ρ), ∞ 1 ρ ρ ∞ [FAGKZ16] M 1, s ∈ [1 + 1/ρ, 1 + ρ), ∞ [1/ρ, s/(1 + ρ)] [ 1 s + ρ 1+ρ , ρ] ∞ [FAGKZ16] M 1, s ≥ 1 + ρ, ∞ 1 1 ∞ [FAGKZ16] LPT M 1, s < ρ, 2 0 ∞ ∞ [FAGKZ16] M 1, s ≥ ρ, ∞ 1 1 ∞ [FAGKZ16] SPT M 1, s < ρ, 2 [ 1 2+ρ , (s-1)ρ +1 (s-1)ρ +1+ρ ] ∞ ∞ [FAGKZ16] M 1, s ≥ ρ, ∞ 1 1 ∞ [FAGKZ16] γ-Burst M 1, [1 + γ/ρ, ρ), 2 1 1 1 [Zav16, KWZ15]
TABLE 1: Detailed metric comparison of online scheduling algorithms for the case of a single machine. The last column provides the references where the results of the corresponding row can be found. Note that by definition, 0-completed-load competitiveness ratio equals to non-competitiveness, as opposed to the other two metrics, where non-competitiveness corresponds to an ∞ competitiveness ratio.

Note that, parameter γ represents the smallest number of c min -tasks that an algorithm running with speedup s can complete, in addition to a c max -task, in an interval of length (γ + 1)c min . What is more interesting, is that in all cases no algorithm can achieve competitiveness neither with respect to pending load nor latency. Then, studying the four popular algorithms mentioned, we prove some positive results, guaranteeing some competitiveness for the different ranges of speedup s and the number of task sizes k. Observe that, with the exception of SPT, none of the algorithms is competitive in any of the three metrics when s < ρ ; algorithm SPT is competitive only in terms of completed load and only when two task sizes are considered. What is more, in terms of latency, only algorithm LIS is competitive when s ≥ ρ. This specific result may not be surprising, since algorithm LIS gives priority to the tasks that have been waiting the longest in the repository. A rather interesting observation though, is that algorithms LPT and SPT become 1-competitive in terms of completed and pending load for s ≥ ρ, whereas LIS and SIS require larger speedup to achieve this. We can say that these results demonstrate some differences between two classes of scheduling policies : the ones giving priority based on the task arrival time (LIS and SIS) and the ones giving priority based on the task size (LPT and SPT). Observe also, that different algorithms scale differently with respect to the speedup ; with the increase of the machine speed the asymptotic competitive performance of each algorithm changes in a different way for each policy. Nonetheless, it is not easy to denote one of the four algorithms as generally better than the rest. The answer would depend on the exact model and objective.

After the limiting results of the widely-used algorithms we propose algorithm γ-Burst, for which we show optimal competitiveness in all three measures for the case of speedup s ∈ [1 + γ/ρ, ρ), however only for two task sizes. The idea of this algorithm is briefly described as follows : If the tasks pending are of the same size, it schedules one of them. Else, if there are at least γ small tasks available, it schedules γ of them and then a large task. Otherwise, it schedules one task of each size interchangeably.

Table 2 is the corresponding table for the results in the case of multiple machines. One can observe right away that we did not consider the latency competitiveness and this is because of its high complexity even in the case of one machine. Nevertheless, we have defined a type of algorithms, called GroupLIS(β), for which we have proven some fundamental properties showing that they avoid redundant task executions. An algorithm is of type GroupLIS if the following three conditions hold : (1) it separates the pending tasks into

M m, 1, 1 1 1 [Zav16] M m, s ≥ ρ, ∞ [1/ρ, 1] [1, ρ] [Zav16] M m, s ≥ 1 + ρ, ∞ 1 1 [Zav16] (m, β)-LIS M m, s ≥ ρ, ∞ [1/ρ, 1] [1, ρ] [FAGKZ15] γm-Burst M m, s ∈ [1 + γ ρ , ρ), 2 1 1 [FAGKZ15] (m, β)-LAF M m, 7/2, k 1 1 [FAGKZ15]
TABLE 2: Detailed metric comparison of the algorithms proposed for the case of multiple machines for different ranges of speedup and number of task sizes. Again, the last column provides the reference where the results of the corresponding row can be found. Note that GroupLIS is a type of parallel algorithms and parameter β is a constant that characterizes the corresponding algorithms.

classes according to their size, (2) it sorts the tasks in each class by their arrival time, and (3) when a class contains at least β • m 2 tasks and machine p decides to schedule a task from that class, then it schedules the (p • βm)th task in the row. We showed that these parallel algorithms actually become optimal for speedup s ≥ 1 + ρ.

Then, trying to find algorithms that would achieve good competitiveness with a lower speedup, we propose (m, β)-LIS, γm-Burst and (m, β)-LAF, each with their limitations. The first two are generalizations of algorithms LIS and γ-Burst presented for the case of one machine. Algorithm (m, β)-LIS actually belongs to the GroupLIS category, so it gives the same performance guarantees. On the other hand, algorithm γm-Burst becomes optimal with a smaller speedup but it only considers two different task sizes. Finally, algorithm (m, β)-LAF considers k finite task sizes and uses an amortization approach to schedule the tasks, but it needs a relatively high speedup, though it could still be preferred depending on the value of ρ. In short, each machine keeps a local variable, named total, where it stores the total load of the completed tasks since their last restart. Each machine then schedules a task from the queue with the largest size such that it is not bigger than parameter total and there are at least β • m 2 tasks pending in that queue.

Discussion

This paper presents only some of the results of my thesis [START_REF] Zavou | Online Scheduling in Fault-prone Systems : Performance Optimization and Energy Efficiency[END_REF] focusing on the most important ones. Some interesting conclusions have been derived with respect to the efficiency of online scheduling algorithms in fault-prone parallel systems in general, as well as for some of the most popular algorithms already used in real life. There are also several questions that remain to be answered, such as the latency study for the multiple machine case. I hope that this work will give the fundamental framework from which to build much further. An interesting extension of this work could be to obtain efficiency bounds as functions of the speedup s used, or try to formalize a single evaluation metric that would encompass the essence of the three we used here. Another, more practical extension of this work, could be to implement real life experiments, for example in data centers, in order to see how the average case results actually scale and whether there is any important aspect in the problem that we currently ignore.

  Completed Load, C Pending Load, P Reference GroupLIS(β)

† This research formed part of my PhD studies in IMDEA Networks and the University of Carlos III of Madrid[START_REF] Zavou | Online Scheduling in Fault-prone Systems : Performance Optimization and Energy Efficiency[END_REF]. I was partially supported by the FPU12/00505 Grant from the Spanish Ministry of Education, Culture and Sports (MECD).