J Yvonnet

D Ryckelynck

P Lorong

F Chinesta

A new extension of the natural element method for non-convex and discontinuous problems: the constrained natural element method (C-NEM)

Keywords: C-NEM, non-convex domains, visibility criterion, constrained Voronoi diagram

In this paper a new extension of the mesh-free natural element method (NEM) is presented. In this approach, coined as constrained natural element method (C-NEM), a visibility criterion is introduced to select natural neighbours in the computation of the shape functions. The computation of these shape functions is based on a modified, constrained Voronoi diagram. With this technique, some difficulties inherent to this method in non-convex domains are avoided and the analysis of problems involving cracks or discontinuities are now easily performed. As the NEM satisfies the Kronecker delta property, the imposition of essential boundary conditions is trivial, unlike other mesh-free methods. The C-NEM technique provides a description of integration cells that allows the use of the stabilized conforming nodal integration (SCNI) scheme instead of Gauss integration to enhance computational efficiency and accuracy. Two numerical examples in elastostatics are reported to evaluate the potential of the proposed technique in highly non-convex geometries, like a crack, through which the solution becomes discontinuous.

INTRODUCTION

In recent years, a new family of computational methods have emerged. The so-called meshless or mesh-free methods have been investigated and used by many researchers for treating a large variety of engineering problems, involving usually large displacements as encountered for example in forming process simulations (free surface or moving boundary problems, moving interfaces, cracks propagation, etc.). In these problems, accurate finite element solutions require significant computational efforts in remeshing steps. In contrast, meshless methods require only nodal data without explicit connectivity between nodes. The first generation of meshless methods, the smooth particle hydrodynamics (SPH) method, originally proposed by Lucy [START_REF] Lucy | A numerical approach to the testing of fusion process[END_REF], introduces kernel approximations and imposes the conservation laws using a collocation technique. Recently, other meshfree methods have been proposed, such as the diffuse element method (DEM) [START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF], the element free Galerkin (EFG) [START_REF] Belytschko | Element-free Galerkin methods[END_REF], the partition of unity method (PUM) [START_REF] Melenk | The partition of unity finite element method: basic theory and applications Computer[END_REF], the reproducing kernel particle method (RKPM) [START_REF] Liu | Reproducing kernel particle methods[END_REF] or the HP-clouds method [START_REF] Duarte | An H-p adaptative method using clouds[END_REF] and they have been successfully used to discretize the Galerkin weak form associated with several problems. Galerkin meshfree methods differ from SPH in two ways: (i) consistency is introduced from kernel corrections, and (ii) discretization is applied to the weak formulation.

Nevertheless, a troublesome task in meshless methods is the imposition of essential boundary conditions because the shape functions related to these methods do not satisfy the Kronecker delta property and they do not vanish on the boundary where the essential conditions are imposed. Special strategies are hence required to enforce the essential boundary conditions, involving additional and costly computational efforts. Among them are the Lagrange multiplier method [START_REF] Belytschko | Element-free Galerkin methods[END_REF], the transformation method [START_REF] Chen | Reproducing kernel particle methods for large deformation analysis of nonlinear structures[END_REF], the approach based on the D'Alembert's principle [START_REF] Günther | Implementation of boundary conditions for meshless methods[END_REF], the introduction of a singular weight function [START_REF] Kaljevic | An improved element free Galerkin formulation[END_REF], the penalty method [START_REF] Zhu | A modified collocation method and a penalty formulation for enforcing the essential boundary conditions in the element free Galerkin methods[END_REF], the mixed transformation method [START_REF] Chen | New boundary conditions treatments in meshfree computation of contact problems[END_REF] and the boundary singular kernel method [START_REF] Chen | New boundary conditions treatments in meshfree computation of contact problems[END_REF]. In order to avoid this additional treatment, a new meshless method has recently been investigated: the natural element method (NEM) [START_REF] Sambridge | Geophysical parameterization and interpolation of irregular data using natural neighbours[END_REF][START_REF] Sukumar | The natural elements method in solid mechanics[END_REF][START_REF] Cueto | On the imposition of essential boundary conditions in natural neighbour Galerkin methods[END_REF]. In the NEM, the test and trial functions are constructed using natural neighbour (Sibson) co-ordinates [START_REF] Sibson | A vector identity for the Dirichlet tesselations[END_REF]. Thus, shape functions satisfy the Kronecker delta property, allowing the imposition of the essential boundary conditions like in the finite element method (FEM). Another important property of the natural neighbour interpolation is that it satisfies linear consistency [START_REF] Sibson | A vector identity for the Dirichlet tesselations[END_REF]. The natural neighbour co-ordinates are based on well-known geometric concepts: the Voronoi diagram and the Delaunay tesselation. As the Voronoi diagram and its dual Delaunay tesselation are only defined for convex domains, issues occur when nonconvex domains are considered. Thus, mutual influences between nodes on the boundary of non-convex domains appear, which introduces difficulties in the imposition of essential boundary conditions.

An extension of the NEM to facilitate its application for any domain (convex or not) and retaining its attractive properties, the constrained natural element method (C-NEM) was proposed in our former work [START_REF] Yvonnet | Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi contraint-Méthode des éléments C-Naturels[END_REF]. The main originality of this approach is the introduction of a visibility criterion [START_REF] Krysl | Element-free Galerkin: convergence of the continuous and discontinuous shape functions[END_REF] in the NEM to select natural neighbours for constructing the interpolation. A modified, so-called constrained Voronoi diagram is introduced for the computation of the shape functions. In addition, the stabilized conforming nodal integration (SCNI) scheme [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF][START_REF] Gonzalez | Numerical integration in natural neighbour Galerkin methods[END_REF] is used instead of Gauss integration to improve both the accuracy and the efficiency of the method.

In this paper, after a brief review of the natural neighbour interpolation and its drawbacks in non-convex domains, the C-NEM is presented. A discussion on its main properties, focusing on its applicability in non-convex domains, is reported. Finally, two numerical examples in elastostatics are presented to illustrate the accuracy of the C-NEM performance in non-convex and highly non-convex domains.

THE NATURAL ELEMENT METHOD

Voronoi diagram and natural neighbours

The notion of natural neighbours and natural neighbours interpolation were introduced by Sibson [START_REF] Sibson | A vector identity for the Dirichlet tesselations[END_REF] for data fitting and smoothing of highly irregular set of nodes. This interpolant is based [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Deuxième mémoire: recherhes sur les parallélloèdres primitifs[END_REF] and later applied in many fields of engineering sciences. In 2D, the Voronoi diagram of an irregular set of nodes divides the plane into a set of regions T i , one for each node n i , such that all points in a particular region are closer to its defining node than to any other node (see Figure 1). The Voronoi diagram is unique, it fills the convex hull of the set of nodes, and can be extended to any dimension. In mathematical terms, the Voronoi diagram is defined as

T i = {x ∈ R n : d(x, x i )<d(x, x j ), ∀j = i}, ∀i (1)
The other closely related geometric construct which is of interest here is the Delaunay tesselation (or triangulation in 2D), which was introduced by Voronoi [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Deuxième mémoire: recherhes sur les parallélloèdres primitifs[END_REF] and extended by Delaunay [START_REF] Delaunay | Sur la sphère vide. A la memoire de Georges Voronoi Izvestiya Akademic[END_REF]. The Delaunay triangles are constructed connecting the nodes whose Voronoi cells have common boundaries. From now on, we use 2D terminology although the discussion is also valid in 3D. The natural neighbours of any node are those defining its neighbouring Voronoi cells, that is, those connected with the node under consideration by the side of a Delaunay triangle.

If the distance between nodes is relatively large in some parts, or its distribution is highly anisotropic, then the set of natural neighbours will reflect these features, but they nevertheless still represents the best set of nearby surrounding nodes. They are natural candidates to define the basis of a local interpolation scheme.

Natural neighbour shape functions

In order to quantify the neighbour relation for any point x introduced into the tesselletion, Sibson [START_REF] Sibson | A vector identity for the Dirichlet tesselations[END_REF] used the concept of second-order Voronoi cells T ij . T ij is defined as the locus of the points x such that its closest node is n i and the second closest node is n j , i.e.

T ij = {x ∈ R n : d(x, x i )<d(x, x j )<d(x, x k ) ∀k = i, ∀k = j, i = j } (2)
Let be a Lebesgue measure (length, area or volume in 1D, 2D or 3D, respectively). Let x be the Lebesgue measure of T x and xi the Lebesgue measure of T xi . The natural neighbour co-ordinates of x with respect to a natural neighbour node n i is defined as the ratio of the Lebesgue measure of second-order cell T xi to the Lebesgue measure of first-order cell T x : If we consider the two-dimensional case shown in Figure 2, the natural neighbour co-ordinate or natural neighbour shape function is given by

i (x) = xi x (3)
i (x) = Area(afghe) Area(abcde) (4)
Braun and Sambridge [START_REF] Sambridge | Geophysical parameterization and interpolation of irregular data using natural neighbours[END_REF] have extended the Sibson theory [START_REF] Sibson | A vector identity for the Dirichlet tesselations[END_REF] to deduce expressions for the derivatives of the interpolation functions, which can be applied in a Galerkin framework for solving partial differential equations, where trial and test functions are interpolated using natural neighbour shape functions. These authors refer to this method as the NEM.

Thus, the interpolation of a vector-valued function u h (x) results

u h (x) = n i=1 i (x)u i ( 5
)
where u i (i ∈ [1, , n]) represents u(x i ) (x i are the co-ordinates of node n i ), n the number of natural neighbour nodes of the point x, and i (x) is the shape function related to node n i and given by Equation (4). 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 neighbours of node n i ). The support of a node n i in a particular nodal distribution is depicted in Figure 3.

Natural element interpolation properties

In this section, the main properties of the natural neighbour shape functions are briefly reviewed. See Reference [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF] for more details.

• The Kronecker delta property: The NEM and finite element shape functions share the following property:

i (x j ) = ij [START_REF] Duarte | An H-p adaptative method using clouds[END_REF] which implies that u h (x i ) = u i . Thus, the imposition of the essential boundary conditions is direct for the concerned nodes, in contrast to most of the other meshless approximations where u h (x i ) = u i . This property and the linear variation of the NEM approximation on convex boundaries allows to easily impose essential boundary conditions. • The partition of unity: By construction we have the following relation:

n i=1 i (x) = 1, ∀x ∈ (7)
This property is equivalent to a zeroth-order consistency, which means that the interpolant can exactly reproduce constant functions. • Linear consistency: The natural neighbour shape functions also satisfy the local co-ordinate property [START_REF] Sibson | A vector identity for the Dirichlet tesselations[END_REF], namely,

x = n i=1 i (x)x i (8)
which combined with Equation ( 7), implies that the natural neighbour interpolant spans the space of linear polynomials (linear completeness). • Linear variation on convex boundaries: The trial function is strictly linear on a Delaunay triangle edge located on the boundary of the convex hull. The proof of this property can be found in Reference [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]. • Smoothness: Natural neighbour 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 (Figure 3) (spheres in 3D) where they are only C 1 , because of the discontinuity on the number of natural neighbours nodes across these boundaries. Hiyoshi and Sughara [START_REF] Hiyoshi | Improving continuity of Voronoi-based interpolation over Delaunay spheres[END_REF] have shown that the Sibson interpolant belong to a more general class of Voronoi-based interpolants, called kth-order standard co-ordinates, proving that the interpolant generated by the kth-order standard co-ordinates have C k continuity on the Delaunay circles (spheres) boundaries. In this context, the Sibsonian and non-Sibsonian (Laplace) co-ordinates [START_REF] Sukumar | Natural neighbour Galerkin methods[END_REF] results to be the standard co-ordinates of order 1 and 0, respectively. cited above are not satisfied anymore. The lack of linear precision on non-convex boundaries is also a consequence of this problem. To illustrate this fact, we consider the non-convex domains depicted in Figure 4, where a point x is introduced. In order to compute the interpolation at that point a new Voronoi cell is created, whereas the Voronoi cells related to the neighbour nodes of point x are modified (see Figure 4). We can see in Figure 4 (left) that the Voronoi cell just created shares edges with the Voronoi cells associated to nodes n 1 , n 2 , n 3 and n 4 that then result natural neighbours of point x. However, in order to impose the essential boundary condition at point x we require both, the shape function associated with node n 3 must vanish at point x and the influence of node n 4 on the solution at point x must be removed, i.e. n 4 cannot be a natural neighbour of point x.

Several techniques have been proposed to overcome this problem lying to the non-convex domains. A coupling between the finite element and natural element method [START_REF] Sukumar | Natural neighbour Galerkin methods[END_REF] or the use of the -shape concept [START_REF] Cueto | Imposing essential boundary conditions in the natural elements method by means of density-scaled alpha-shapes[END_REF] can be cited among many other possibilities.

However, some difficulties persist in the case of highly non-convex domains, as encountered in the crack analysis for example. In the context of the -NEM, Cueto et al. [START_REF] Cueto | On the imposition of essential boundary conditions in natural neighbour Galerkin methods[END_REF] propose to describe the geometry of the domain in -shape sense, which requires a control of the nodal density. Furthermore, the same authors have shown that it is not possible to describe by this way sharp corner or cracks. On the other hand, Sukumar [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF] propose the introduction of an explicit description of the domain boundary by means of a planar straight graph (PSG). We propose in the next section a simple extension of the NEM to circumvent the just referred difficulties.

THE CONSTRAINED NATURAL ELEMENT METHOD (C-NEM)

In this section, a new extension of the NEM is proposed, which makes use of a visibility criterion, introduced to restrict natural neighbours. A constrained Voronoi diagram results from the application of that criterion, which will be used to define the interpolation functions, keeping, as later discussed, the main properties of the NEM interpolation.

The visibility criterion

In our opinion, the simplest way to solve the difficulties discussed in the previous section lies in the introduction of a criterion to remove parasite influences between nodes on non-convex boundaries. The criterion must interact with a description of the boundary in a CAD sense.

Organ et al. have introduced in Reference [START_REF] Organ | Continuous meshless approximations for nonconvex bodies by diffraction and transparency[END_REF] a visibility criterion to define the domain of influence of a node near a crack tip and other non-convex domains. In their approach, the point x is in the domain of influence of a node n i if x is within the region where the shape function associated with n i is non-zero ( i >0) and it is visible from node n i when the boundaries are assumed opaque. Thus, we can define

i = {x | x ∈ , i (x)>0} V IS i = {x | S x → x i ∩ = ∅} (9)
where S x→x i represents the straight line connecting x and n i , i.e.

S x→x i = x + (x i -x), ∈ [0, . . . , 1] (10)
where x i is the vector defining the position of node n i and is the domain boundary (including interior parts). Thus, the domain of influence of a node n i , * , is defined by

* i = VIS i ∩ i (11)

C-NEM interpolation

If the visibility criterion is introduced in the NEM, the natural neighours become constrained natural neighbours (C-n-n). The set of natural neighbours will be restricted by applying the visibility criterion introduced in the previous section. Thus, the new functional approximation can be written as

u h (x) = V i=1 C i (x)u i ( 12
)
where V is the number of natural neighbours visible from point x and C i is the constrained natural neighbour shape function. The computation of the C-n-n shape function is similar to natural neighbour shape function, when one proceed on the constrained Voronoi diagram introduced later.

Support of constrained natural neighbour shape functions

In the C-NEM, the support of C i consists of the union of the n circles (spheres in 3D), passing through the vertices of the Delaunay triangles (tetrahedra in 3D) defined by the node under consideration n i , and the V natural neighbours visible from node n i . Furthermore, as stated by Edelsbrunner [START_REF] Edelsbrunner | Geometry and Topology for Mesh Generation[END_REF], these circles can extend across the boundaries, but cannot contain nonvisible neighbours. In Figure 5, the intersection between the support of a constrained natural neighbour shape function C i and the boundary is depicted.

Γ Ω Γ Ω Figure 5
. Supports of the natural neighbour and constrained natural neighbour shape functions.

Sibson [START_REF] Sibson | A vector identity for the Dirichlet tesselations[END_REF] proved that the partition of unity and the linear consistency are satisfied in the context of natural neighbour interpolants, if, and only if, the point where the interpolation is computed lies inside the convex hull of their natural neighbours. It can be noted that in the context of the constrained natural neighbours interpolation, any point in the domain lies inside the convex hull of their constrained (visible) natural neighbours. In consequence, the partition of unity and the linear consistency are conserved in the C-NEM.

Constrained Voronoi diagram

In this section, the constrained Voronoi diagram (CVD) is introduced to provide a suitable background for the computation of natural neighbour shape function over any kind of geometry (convex or not). The constrained Voronoi diagram (or extended) Voronoi diagram, is introduced by Seidel in Reference [START_REF] Seidel | Constrained Delaunay triangulations and Voronoi diagrams with obstacles[END_REF]. In our work, we prefer the term constrained instead of extended to avoid confusions with the extended finite element terminology [START_REF] Melenk | The partition of unity finite element method: basic theory and applications Computer[END_REF]. We consider the definition of the CVD provided in Reference [START_REF] Edelsbrunner | Geometry and Topology for Mesh Generation[END_REF]: Imagine R 2 is a sheet of paper 0 , with the points of the set of nodes S and the line segments defining the boundary in a set L, drawn on it. For each l i ∈ L, we cut 0 open along l i and glue another sheet i , which also cut open along l i . The gluing is done around l i such every traveller who cross l i switches from 0 to i and vice versa. A schematic view of the particular gluing necessary to achieve that effect is illustrated in Figure 6(a). We know what it means for two points on the primary sheet to be visible from each other. For other pairs we need a more general definition. For i = 0, points x 0 ∈ 0 and y i ∈ i are visible if xy crosses l i , and l i is the first constraining line segment crossed if we traverse xy in the direction from x to y.

In Figure 6(b), the intersection between the CVD and the domain closure is depicted. The resulting diagram is composed with cells T C i , one for each node n i , such as any point x inside T C i is closer to n i than to any other node n j visible from x. We call it the constrained Voronoi cells, which are defined formally by

T C i = {x ∈ R n : d(x, x i )<d(x, x j ), ∀j = i, S x→n i ∩ = ∅, S x→n j ∩ = ∅} (13)
where is the domain boundary, composed with segments l i ∈ L. A generalization of the constrained Delaunay triangulation to 3D does not exist without adding new nodes, as shown in Reference [START_REF] Schönhardt | Uber die zerlegung von dreieckspolyedern in tetraeder[END_REF]. Nevertheless, some techniques to construct 3D constrained Delaunay tesselations are provided in References [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 addition of Steiner points. The constrained Voronoi diagram has been introduced for two reasons:

(i) Once such diagram is constructed, classical algorithms for the computation of the shape functions, Watson-Bowyer [START_REF] Sambridge | Geophysical parameterization and interpolation of irregular data using natural neighbours[END_REF] or Lasserre [START_REF] Sambridge | Computational methods for natural neighbour interpolation in two and three dimensions[END_REF] can be applied directly because connections between non-visible natural neighbours in non-convex domains have been removed. (ii) Recently a new integration method, the SCNI [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] has been introduced in the meshless methods framework for improving both the accuracy and the efficiency of the numerical integration of the weak formulation. This technique, as proved later, is specially appropriate in the context of NEM, because the integration cells match the constrained Voronoi cells.

Computation of the C-n-n shape functions

Once the constrained Voronoi diagram is computed, selection of constrained natural neighbours is a direct task because this diagram only contains connections between mutually visible nodes. In Figure 7(a), an example of the computation of a natural shape function making use of the standard Voronoi diagram is presented. We can notice that the n 1 , n 2 n 3 and n 4 nodal values will contribute to the interpolation in x (located on the non-convex domain boundary). We can observe that the computed solution in x is not good because its value does not result from a linear interpolation between the nodal values n 1 and n 3 . Thus, the imposition of a essential boundary condition is disturbed by the influence of internal (n 2 ) and external and consequently non-visible (n 4 ) nodes.

In Figure 7(b), the same computation is done making use of the constrained Voronoi diagram. In such a diagram, as previously described, any Voronoi cell facet related to two nodes which are not mutually visible are removed. In this case, the interpolated value in x will depend only on nodal values in n 1 , n 2 and n 3 . Moreover, taking into account the fact that the second-order Voronoi cells T 1x and T 3x have infinite areas, the contribution of the internal node n 2 vanishes (the area of T 2x is finite).

We can summarize these facts in the following manner: the second-order Voronoi cells related to nodes located on the domain boundary (convex or not), have an infinite area because the influence between non-visible nodes has been removed during the construction of the constrained Voronoi diagram. In that case, influences of interior nodes also vanish. Thus, essential boundary conditions can be easily enforced on any domain boundary (convex or not), where the approximation is in fact linear (Figure 8).

Remark 1

To compute a general constrained Voronoi diagram only a description of the boundary in a CAD sense is required. We can also introduce inclusions or discontinuities from a list of edges (facets in 3D) defining their boundaries. Crack analysis with eventual auto-contact can be then also treated as well as moving boundaries, interfaces or discontinuities when its kinematics is known.

Remark 2

New nodes can be added (or removed) in the domain or/and on the domain boundary in order to improve the solution accuracy or to achieve an accurate representation of moving boundaries. In the nodes just introduced a C-NEM interpolation of internal variables must be carried out. In the context of the -NEM, it is required that the nodal cloud of points must reproduce the geometry of the domain in an -shape sense, to remove spurious influences between nodes and to ensure strict linear interpolation along the boundaries. This is not always possible, as shown in Reference [START_REF] Cueto | On the imposition of essential boundary conditions in natural neighbour Galerkin methods[END_REF].

NUMERICAL EXAMPLES

Elastostatic formulation

We consider the 2D small displacement elastostatics problem, which is described by the equilibrium equation:

∇ • + b = 0 in ( 14
)
where ∈ R 2 is the material domain, ∇ is the gradient operator, is the Cauchy stress tensor and b is a body force term. The constitutive relation is given by = C [START_REF] Sibson | A vector identity for the Dirichlet tesselations[END_REF] where and are the vector form of the stress and linearized strain tensor (symmetric part of the gradient of displacements tensor), respectively, which are written using the conventional notations:

=    11 22 12    , =    11 22 2 12    ( 16
)
and C is the elastic tensor matrix form.

The essential and natural boundary conditions are

u = ū on u n = t on t ( 17
)
where = u ∪ t is the boundary of , n is the unit outward vector defined on , and ū and t the prescribed displacements and tractions, respectively.

The variational formulation (principle of virtual work) associated with the elastostatic problem results:

Find u ∈ H 1 ( ) kinematically admissible

(u = ū on u ) such that * • d = v * • b d + t v * • t d , ∀v * ∈ H 0 1 ( ) ( 18
)
where H 1 ( ) and H 0 1 ( ) are the usual Sobolev functional spaces. On substituting the trial and test functions (both approximated in the C-NEM framework) in the above equation and using the arbitrariness of the field v * , the following linear system of equations is obtained after numerical integration:

Kd = f ext ( 19
)
where d is the vector containing the nodal displacements, and the matrix K is given by

K = B t CB d ( 20
)
where C is the usual matrix form of the elastic tensor

f ext = t N t t d + N t b d (21)
N is the matrix containing the shape functions:

N = 1 0 2 0 . . . N 0 0 1 0 2 . . . 0 N ( 22
)
B is the matrix containing the shape functions derivatives:

B =     1,x 0 2,x 0 . . . N,x 0 0 1,y 0 2,y . . . 0 N,y 1,y 1,x 2,y 2,x . . . N,y N,x     ( 23
)
The constitutive matrix C in isotropic linear elasticity results for plane strain:

C = E (1 -2 )(1 + )     1 - 0 1 - 0 0 0 1-2 2     (24)
and in the case of plane stress:

C = E (1 -2 )     1 0 1 0 0 0 1- 2     (25)

Stabilized numerical integration

Recently, new numerical integration procedures have been proposed in the framework of meshless simulations. Domain integration using Gauss quadrature introduces significant numerical error when integration cells do not match shape function supports. Another source of numerical errors is the non-polynomial form of the shape functions that makes non-optimal the Gauss integration [START_REF] Dolbow | Numerical integration of the Galerkin weak forme in meshfree methods[END_REF]. On the other hand, direct nodal integration, using the nodes as integration points, leads to numerical instabilities.

In this work, we use the integration procedure proposed by Chen et al. [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF]: the stabilized conforming nodal integration (SCNI). In this approach, a strain smoothing stabilization is performed to stabilize the nodal integration.

The SCNI procedure is based on the assumed strain method, in which a modified gradient is introduced at the integration point (node) [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF]:

∇u h j (x k ) = 1 A k k ∇u h j (x k ) d (26)
where x k are the co-ordinates of node n k . Thus, the modified strain vector is given by

˜ h (x k ) = 1 A k k h (x) d = 1 A k k      *u h *x *v h *y *u h *y + *v h *x      d (27)
Now, by applying the divergence theorem, the following equation is deduced:

˜ h (x k ) = 1 A k k     u h n x v h n y u h n y + v h n x     d ( 28
)
where k is a representative domain related to node n k and k and A k are its boundary and area, respectively, as depicted in Figure 9.

Introducing now the natural element approximation, i.e. u h and v h according to Equation ( 12) we obtain which is expressed explicitly by

˜ h (x k ) = Bk d (29)
n k Ω k Γ k (a) (b)
˜ h (x k ) =     bx1 (x k ) 0 bx2 (x k ) 0 . . . bxN (x k ) 0 0 by1 (x k ) 0 by2 (x k ) . . . 0 byN (x k ) by1 (x k ) bx1 (x k ) by2 (x k ) bx2 (x k ) . . . byN (x k ) bxN (x k )                     u 1 v 1 u 2 v 2 . . . u N v N                 = Bk d (30)
where bxj (x k ) and byj (x k ) are defined by bxj

(x k ) = 1 A k k j (x)n x (x) d (31) byj (x k ) = 1 A k k j (x)n y (x) d (32)
Obviously, a lot of components in the matrix Bk are zero due to the compact support of the shape functions. Chen et al. [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] propose to use the intersection between the Voronoi diagram and the domain closure to define the integration cells k . Introducing the C-NEM approximation, the smoothed strain and the nodal integration, the global stiffness matrix is obtained by assembling the contribution of each node n k :

R int R ext R 1 Ω 2 Ω 1 Γ 1 Γ 2 (a) (b)
K = k Kk = k A k Bt k C Bk (33)
The stress in each cell is then deduced from

˜ h k = C Bk d (34)
The application of the divergence theorem avoid the computation of the shape functions derivatives, and only the evaluation of the shape functions on the boundary of the Voronoi cells is required. The stabilized nodal integration has been successfully applied by Gonzalez et al. [START_REF] Gonzalez | Numerical integration in natural neighbour Galerkin methods[END_REF] in the context of the -NEM, where a significant improvement in the accuracy was noticed.

This new integration technique can be naturally applied in the C-NEM where the nodal integration domains coincide with the intersection between the constrained Voronoi cells and the domain closure. It can be seen in Figure 9 that the intersection between the standard Voronoi diagram and the domain closure does not define accurately the integration domain.

A first example: a two-materials cylinder under internal pressure

We consider firstly the two-layers hollow cylinder under internal pressure whose geometry is depicted in Figure 10 (symmetry conditions has been taken into account. The main interest of this elastostatic simulation is the analysis of the C-NEM capabilities to treat non-convex geometries, handling eventual interface discontinuities. Four clouds of nodes have been used, containing 147, 315, 687 and 1863 nodes. Both regular and irregular distributed nodes were tested.

The exact displacements solution of a hollow cylinder under internal pressure results

u r = ar + b r (35)
In our case we can write the displacements in each domain in the form

u r = a 1 r + b 1 r R int r R 1 ( 36
)
and

u r = a 2 r + b 2 r R 1 r R ext ( 37
)
from which the strains and stresses can be deduced

1 rr = a 1 - b 1 r 2 2 rr = a 2 - b 2 r 2 1 = a 1 + b 1 r 2 2 = a 2 + b 2 r 2 1 r = 0 2 r = 0 (38) 1 rr = a 1 ( 1 + 2 1 ) - 2b 1 1 r 2 2 rr = a 2 ( 2 + 2 2 ) - 2b 2 2 r 2 1 = a 1 ( 1 + 2 1 ) + 2b 1 1 r 2 2 = a 2 ( 2 + 2 2 ) + 2b 2 2 r 2 1 r = 0 2 r = 0 (39)
where i and i , i = 1, 2, denote the Lame's coefficients of both materials, related to elasticity modulus and Poisson's coefficient by

= E 2(1 + ) and = E (1 -2 )(1 + ) (40)
The parameters a i and b i , i = 1, 2, can be computed by applying the following boundary and transmission conditions assuming plane stress:

-rr (R int ) = P rr (R ext ) = 0 1 rr (R 1 ) = 2 rr (R 1 ) u 1 r (R 1 ) = u 2 r (R 1 ) (41)
where P represents the internal pressure. Numerical results of displacements and stress are compared to the analytical solution for a cloud of 315 regularly distributed nodes in Figures 11 and12 (R int = 0.5 mm, R 1 = 0.7 mm, R ext = 0.9 mm, 1 = 0.3, E 1 = 210 000 MPa, 2 = 0.2, E 2 = 105 000 MPa, P = 1 MPa). The C-NEM is in good agreement with the analytical solution, and the discontinuity across the interface is accurately reproduced. Figure 13 and Table I group the errors using the energy norm for the different clouds of nodes and both regular and irregular nodal distributions, where ˜ h is given by Equation [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]. In Figure 13, h denotes the radial distance between two nodes.

A second example: a mode-I crack analysis

A crack analysis using the NEM can be found in Reference [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF], where due to the symmetry of the problem the computation domain is actually convex. It shows that the NEM solution convergence is faster than the finite element one. In the present example, a mode-I cracked plate is considered. The geometry of the model is depicted in Figure 14(a). The analytical solution in the neighbourhood of the crack tip is given by [START_REF] Perrone | Fracture Mechanics[END_REF] 11 =

K I √ 2 r cos 2 1 -sin 2 sin 3 2 (42) 22 = K I √ 2 r cos 2 1 + sin 2 sin 3 2 (43) 12 = K I √ 2 r sin 2 cos 2 cos 3 2 (44)
In order to investigate the convergence of the numerical solution in the whole domain we enforce the boundary conditions F d = ex n on the boundary of the problem domain.

In this example, L = 1 mm, a = L/2 mm and K I = 1. In order to illustrate the ability of the C-NEM to handle the crack behaviour, no symmetry was considered. The convergence of the C-NEM solution is depicted in Figure 15. The parameter h denotes the distance between two nodes in a regular set of nodes. We can notice that the convergence rate is 0.42, which is approximatively the expected convergence rate (0.5) obtained with triangular linear finite elements.

A third example: a plate containing an internal crack

In this example, we consider the problem defined in Figure 16 where 2a = 0.2 mm, 2W = 2 mm and ∞ 22 = 1 MPa. Plane stress conditions are assumed with an elastic modulus of E = 1 MPa and a Poisson coefficient of = 0.3. Firstly, we consider the case of = 0 whose exact solution in an infinite plate was obtained by Muskelishvili [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]. For = 0 and = 0 in

( = 0, = 0, r) = ∞ 22 a + r √ r(2a + r) 11 ( = 0, = 0, r) = ∞ 22 -22 (45)
When the crack size is much smaller than the domain size, Muskelishvili's solution can be taken as a reference solution (Figure 17).

A cloud of 760 nodes is considered in the C-NEM simulation (400 nodes are uniformly distributed, other 20 nodes are located on each crack side and the remaining 320 nodes are concentrated in the crack tip zones) as depicted in Figure 16 without assuming any symmetry condition. Numerical results are reported in Figure 18. Now, we consider the case of an inclined crack. The stress intensity factor (SIF) K I has been computed using the contour form of the interaction integral method [START_REF] Yau | A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity[END_REF], the different values of interest being easily interpolated anywhere in the domain. The size of the square domain used in the computation is 0.6a × 0.6a (see Figure 6). Numerical results concerning the intensity factor K I for different values of are compared in Figure 19 with the exact ones given in Reference [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]:

K I = ∞ 22 √ a cos 2 (46)
Accuracy can be improved by increasing the nodal density in the crack tip neighbourhood. Another possibility is to extrinsically enhance the trial function by adding terms reproducing the local √ r behaviour. As the shape functions C i (x) define a partition of unity, the enrichment of the C-NEM trial function in crack problems can be achieved in the partition of unity framework [START_REF] Melenk | The partition of unity finite element method: basic theory and applications Computer[END_REF]. Let 0 ⊂ be a sub-domain in a region around the crack tip where the crack tip singularity dominates. Then the enriched function for two-dimensional crack problems could follow the form suggested by Sukumar in Reference [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]:

u h (x) = n i=1 c i (x)u i + √ r m j =1 c j (x)a j ( 47
)
where a j are additional nodal coefficients that are associated with the m nodes in 0 .

CONCLUSION

This paper proposes a new extension of the natural element method, the constrained natural element method (C-NEM) applicable to any geometry (convex or not) involving eventual discontinuities (fixed or moving). For this purpose, a visibility criterion has been introduced to avoid undesirable nodal influences in non-convex domains. The constrained or extended Voronoi diagram which is known as the strict dual to the constrained Delaunay tesselation, is introduced as the basis for computing the interpolation functions and it is also applied for defining the integration domains used in the stabilized nodal integration (SCNI) [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF][START_REF] Gonzalez | Numerical integration in natural neighbour Galerkin methods[END_REF]. The main properties of the C-NEM can be summarized as follows:

(i) Direct enforcement of essential boundary conditions without any special treatment.

(ii) The ability to reproduce discontinuities across interfaces.

(iii) Continuity inside the domain, which is not perturbed by the introduction of the visibility criterion. (iv) Linear precision on any boundary (convex or not), without parameter dependance. (v) A significant improvement of both efficiency and accuracy in the numerical integration by using the SCNI procedure, as demonstrated in Reference [START_REF] Gonzalez | Numerical integration in natural neighbour Galerkin methods[END_REF]. (vi) Linear consistency and partition of unity property at point of the domain. The second property allows functional enrichment.

Topics related to enrichment and the natural extension of the C-NEM for treating 3D models, constitute some of our current work in progress.

Figure 1 .

 1 Figure 1. Voronoi diagram of a set of nodes in the plane.

Figure 2 .

 2 Figure 2. Construction of natural neighbour shape functions.

Figure 3 .

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

Figure 4 .

 4 Figure 4. Influences from nodes lying on non-convex boundaries.

Figure 6 .

 6 Figure 6. Interpretation of the constrained Voronoi diagram: (a) schematic view of the CVD; and (b) intersection between CVD and the domain closure.

Figure 7 .

 7 Figure 7. Computation of natural neighbour shape functions over a non-convex domain: (a) over standard Voronoi diagram; and (b) over constrained Voronoi diagram.

Figure 8 .

 8 Figure 8. View showing the linearity of the C-NEM shape function associated with node n i along a discontinuity or a non-convex boundary.

Figure 9 .

 9 Figure 9. Integration domains for the SCNI procedure: (a) intersection between the Voronoi diagram and the domain closure; and (b) constrained and bounded Voronoi diagram.

Figure 10 .

 10 Figure 10. (a) Model geometry of bimaterial cylinder; and (b) constrained Voronoi diagram of bimaterial cylinder with irregular spaced nodes.

θFigure 11 .

 11 Figure 11. Displacements across the interface of the bimaterial cylinder.

Figure 12 .

 12 Figure 12. Stresses across the interface of the bimaterial cylinder.

Figure 13 .

 13 Figure 13. Rate of convergence of the C-NEM for the bimaterial cylinder problem using a regular nodal distribution.

Figure 14 .Figure 15 .

 1415 Figure 14. (a) Geometry for the mode-I crack problem; and (b) exact boundary conditions.

Figure 16 .

 16 Figure 16. Square plate with an inclined crack: (a) geometry and loads; and (b) distribution of nodes.

Figure 17 .

 17 Figure 17. Zoom of the constrained Voronoi diagram in the neighbourhood of the crack (detail).

Figure 16 ,

 16 Figure 16, this solution results 22 ( = 0, = 0, r) = ∞

StressesFigure 18 .Figure 19 .

 1819 Figure 18. Radial stresses ahead of the crack tip for = 0.

Table I .

 I Error in the energy norm for the bicylinder problem.

		˜ h -ex E	˜ h -ex E
	Number of nodes	Regular nodal distribution	Irregular nodal distribution
	147	1.13 × 10 -4	1.21 × 10 -4
	315	7.41 × 10 -5	8.10 × 10 -5
	687	4.95 × 10 -5	5.48 × 10 -5
	1863	3.06 × 10 -5	3.46 × 10 -5

2.4. Issues on non-convex domainsSince the natural neighbour interpolant is constructed on the basis of the Voronoi diagram and Delaunay tessellation, issues occur when non-convex domains are considered. The Voronoi diagram is constructed filling the convex hull of the cloud of nodes, without taking into account the actual domain boundary.Unfortunately, on non-convex domains, some nodes lying on non-convex boundaries can share a Voronoi cell facet and thus they become natural neighbours. In this situation, the properties

The discontinuity is explicitly described by nodes. Owing to the salient shape of the support, which passes through the visible neighbour nodes, the shape functions are strictly continuous near a crack tip, without requiring the use of particular techniques to ensure the continuity[START_REF] Organ | Continuous meshless approximations for nonconvex bodies by diffraction and transparency[END_REF].


D R J Chillingworth

A B Nordmark
email: nordmark@mech.kth.se

Periodic Orbits Close to Grazing for an Impact Oscillator

We show how the geometric impact surface approach to the dynamics of an impact oscillator provides an immediate visualization of the criteria that determine the existence of an impacting periodic orbit close to grazing. We recover the criteria set out earlier by A. Nordmark and indicate how the geometric setting and singularity geometry may be exploited to yield appropriate criteria in degenerate situations where the Nordmark criteria would not apply.

The Impact Oscillator

, for examples and further references.

The simplest model for such a system is the one degree of freedom impact oscillator of the form

where f : R 2 → R and g : R → R are smooth functions, g is T -periodic with T > 0, and where x is restricted to the region x ≥ c. To describe the local

consequence of impact with the obstacle at x = c a restitution law is assumed, which usually takes the form of replacing ẋ < 0 by -r ẋ, r > 0 at impact. Impact oscillators can be seen as special cases of N -degree of freedom (autonomous) systems with codimension-1 obstacles, or of general hybrid systems involving piecewise-smooth vector fields in R n together with maps defined on the codimension-1 boundaries of the regions of smoothness: see [START_REF] Di Bernardo | Piecewise-Smooth Dynamical Systems: Theory and Applications[END_REF] for further discussion.

An impact oscillator (2.1) exhibits certain bifurcation phenomena that arise specifically from the impacts, in addition to those seen in smooth (non-impacting) systems. The most natural example to study first is that in which a hyperbolic (or, at least, isolated and persistent) periodic orbit of period T lying in the region x > c collides with the obstacle at x = c as a parameter μ in the system increases through a critical value μ 0 . For μ = μ 0 the orbit satisfies ẋ = 0 where x = c, that is it grazes the obstacle. It is not a priori clear whether such a periodic orbit will persist (with impacts) for μ > μ 0 or how the result may depend on the nature of the Poincaré map of the grazing orbit.

There have been many studies, both theoretical and numerical, of the dynamical behaviour of an impact oscillator close to grazing: see for example [2, 4, 6-8, 10-12, 15, 16]. A thorough investigation of the persistence or otherwise of a grazing T -periodic orbit was carried out by Nordmark in [START_REF] Nordmark | Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators[END_REF] in the wider context of an n-degree of freedom (autonomous) system with a codimension-1 obstacle, and as an illustration in that paper the general results are applied to the 1-degree of freedom system (2.1). A detailed analysis of criteria for the persistence of orbits of period nT for n = 1, 2, 3 . . . is carried out, leading to interesting algebra related to the eigenvalues of the Poincaré map. However, simple criteria in terms of the Poincaré map itself are not immediately visible.

The purpose of this paper is, first, to extract the essential data from the general formalism of [START_REF] Nordmark | Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators[END_REF] in order to write down what is needed for the impact oscillator (2.1) in the cases n = 1 and n = 2. Having done this, we then take a different approach and use the impact surface description formulated in [START_REF] Chillingworth | Discontinuity geometry for an impact oscillator[END_REF] in order to reveal in this setting the extremely simple geometric criteria that underlie the phenomenon of grazing bifurcation. Finally we show how this geometry, when converted into algebra, yields the Nordmark criteria.

Nordmark's Criteria

We apply the general formalism from [START_REF] Nordmark | Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators[END_REF] to this setting, and in particular to the cases designated there as (n) where n = 1 or n = 2, that is cases of impacting orbits close to grazing that have a single impact during a period of T or of 2T . The aim is to derive simple criteria for the existence of such periodic orbits as a system with a grazing orbit P of period T is perturbed by a parameter μ. We suppose without loss of generality that the grazing occurs when μ = 0.

Single Impact T -Periodic Orbits

Following the notation from [START_REF] Nordmark | Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators[END_REF], let

A = α γ β δ
denote the Jacobian matrix of the time-T map in the (x, ẋ)-plane for the system (2.1) (in the absence of the obstacle) at the point p = (c, 0) when μ = 0.

The T -periodicity of P means that p is a fixed point of the map. It is convenient to introduce also the covector C = [1, 0] and the vector B = [0, 1] t where t denotes transpose. 1First assume that

Δ -:= det(I -A) = 0
so that P (although not necessarily hyperbolic) persists as a T -periodic orbit P μ in the absence of the obstacle for small |μ| > 0. A central role in [START_REF] Nordmark | Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators[END_REF] is played by the quantity

e : = C(I -A) -1 (κ, λ) t (2.2) = Δ -1 -(1 -δ)κ + γλ (2.3)
where (κ, λ) = ( ∂x ∂μ , ∂ ẋ ∂μ ) evaluated at p when μ = 0. To first order in μ, the location of a T -periodic (n = 1) single-impacting orbit (when it exists) relative to the grazing point p is given by μX ∈ R 2 where ([13, eq. (51)])

C(I -A)X = 0 , CX = -1 (2.4)
and from [13, eq. ( 34)] the condition for such an orbit to exist is

y 1 eμ > 0 (2.5)
where from [13, eq. (32)] we have (up to a positive multiple)

y 1 = B t (I -A)X. (2.6)
From the second equation of (2.4) we see X has the form X = (-1, x 1 ) t , and then combining (2.4) and (2.6) we may write

(I -A) -1 0 x 1 1 = 0 -γ y 1 1 -δ (2.7)
which on taking determinants gives

γ y 1 = -Δ -. (2.8)
Assuming γ = 0 this gives

y 1 e = -γ -1 (1 -δ)κ + γλ
and so the criterion (2.5) gives the following basic result:

Proposition 2.1 (Nordmark). Suppose Δ -= 0 ,γ = 0 and also

ψ := (1 -δ)κ + γλ = 0.
For |μ| sufficiently small, a single-impact orbit of period T exists close to the grazing orbit if and only if the parameter μ has the opposite sign from γψ.

In particular this implies that as μ passes through zero an impacting period-T orbit must exist on one side or the other. If it exists for μ > 0 then the impacting orbit is a continuation of the non-impacting orbit, while for the case μ < 0 the free and the impacting orbit exist simultaneously and mutually annihilate at μ = 0. The latter phenomenon is often called a nonsmooth fold: see [START_REF] Di Bernardo | Piecewise-Smooth Dynamical Systems: Theory and Applications[END_REF].

Single Impact 2T -Periodic Orbits

For a 2T -periodic orbit (n = 2) with single impact the location X is now by analogy with (2.4) given by

C(I -A 2 )X = 0 , CX = -1 (2.9)
while the conditions for such an orbit to exist are b 2,1 e μ > 0 (2.10)

y 2 e μ > 0 (2.11)
where b 2,1 = CAX + 1 (2.12)

y 2 = B t (I -A 2 )X. (2.13) Let α 2 , γ 2 denote the α, γ-entries in A 2 , that is α 2 = α 2 + βγ = αη -det A (2.14) γ 2 = γ η (2.15)
where η = tr A is the trace of A. By analogy with (2.7) we have on taking determinants

ηγ y 2 = -det(I -A 2 ) = -Δ -Δ + (2.16)
with Δ + = det(I + A), and so

γy 2 e = -η -1 Δ + (1 -δ)κ + γλ provided η = 0. If (2.11
) holds then the condition (2.10) can simply be expressed by saying that b 2,1 and y 2 have the same sign. Now if X = (-1, x 2 ) t satisfies (2.9) we find

b 2,1 = -α + γx 2 + 1 (2.17) = -α + γγ -1 2 (α 2 -1) + 1 (2.18) = -η -1 Δ - (2.19)
using (2.14) and (2.15). From (2.16) the same-sign condition on b 2,1 and y 2 is therefore precisely the condition that γ Δ + > 0. To summarize:

Proposition 2.2 (Nordmark).

Suppose Δ -= 0, ηψ = 0 and also γΔ + = 0.

For |μ| sufficiently small, a single-impact orbit of period 2T exists close to the grazing orbit if and only if γΔ + > 0 and μ has the opposite sign from ηψ.

The Impact Surface Approach

We now turn to address the same questions, but from the point of view of the geometry of the impact surface as set out in [START_REF] Chillingworth | Discontinuity geometry for an impact oscillator[END_REF]. We briefly describe the formalism, and refer to [START_REF] Chillingworth | Discontinuity geometry for an impact oscillator[END_REF] for further details.

The impact surface V c is defined as

V c = {(τ, v, t) ∈ R 3 : x c (τ, v, t) = c} (2.20)
where x(c, v, τ ; u) denotes the solution of (2.1) with initial data (x, ẋ) = (c, v) when u = τ , and we write

x c (τ, v, t) = x(c, v, τ ; τ + t).
Clearly V c contains the t = 0 plane Π, and in [START_REF] Chillingworth | Discontinuity geometry for an impact oscillator[END_REF] it is shown that V c = Π ∪ V c where V c is a smooth 2-manifold (except generically for finitely many values of c at which V c undergoes a Morse perestroika) that intersects Π transversely along the τ -axis. The impacting system can then be modelled by a discontinuous discrete dynamical system on Π defined as the composition

G c := R • φ c • F c : Π → Π where F c : Π → V c is the first hit map (τ, v) → (τ, v, t 1 ) where t 1 is the smallest positive value of t for which (τ, v, t) ∈ V c φ c : V c → Π is the reset map (τ, v, t) → (τ + t, ẋ(τ + t)) R : Π → Π is the restitution map (such as (τ, v) → (τ, -rv)).
In other words, we reconstruct the dynamics by taking each impact at x = c as a new set of initial data and then proceeding (after applying restitution) to the next impact.

Remark 2.1. In general the map F c may fail to be defined for either of the following two reasons: for given (τ, v) ∈ Π the subset {t : (τ, v, t) ∈ V c } of R may be empty, or (when v = 0) may accumulate at zero. Neither of these arises in the present discussion, however. It is also the case that infinitely many impacts can occur in a finite time interval (the chatter phenomenon), but again for current purposes this does not concern us.

Singularity Geometry

The singularity structure of the projection map

π c = π|V c : (τ, v, t) → (τ, v)
and of the reset map φ c are crucial in construction of the dynamical map G c : Π → Π. For (τ, 0, t) ∈ V c let a = a(τ ) := ẍc (τ, 0, 0) denote the 'initial' acceleration. Then from [START_REF] Chillingworth | Discontinuity geometry for an impact oscillator[END_REF] we have Proposition 2.3.

The singular set S(π c

) consists of the points on V c where ẋc = 0.

If ẍc = 0 then the singularity is a fold singularity. 2. The singular set S(φ c ) consists of the points on V c where v = 0.

If a = 0 then the singularity is a fold singularity.

The first part of the proposition is elementary, from the geometry of V c itself, although the second part is less obvious. Note that π c S(π c ) is the apparent outline (or contour) of V c as viewed along the t-axis; in [START_REF] Chillingworth | Discontinuity geometry for an impact oscillator[END_REF] we call S(π c ) the horizon curve H. Consider now a grazing and otherwise non-impacting periodic orbit P with period T . In terms of the impact surface the orbit P is represented by a point p 0 = (τ 0 , 0, 0) on the τ -axis such that the line l(τ 0 ) through p 0 parallel to the t-axis intersects V c precisely at the points p n = (τ 0 , 0, nT ), n ∈ Z and is tangent to V c there. Suppose the grazing is quadratic, so that a > 0 and ẍc > 0 at all points p n . From Proposition 2.3 we obtain Proposition 2.4. Each of the maps π c : V c → Π and φ c : V c → Π has a fold singularity at each of the points p n , n = 0. See Fig. 2.1 for an illustration of this geometry for n = 1.

The T -periodicity of g means that the impact surface V c is invariant under translation τ → τ + T , so working modulo T on the τ -axis we can regard φ c as taking every p n ∈ V c to p 0 ∈ Π, as does π c . Therefore at each p n , n = 0, we have two fold singularities

π c , φ c : (V c , p n ) → (Π, p 0 ).
We now use this geometry to derive a criterion for the existence of impacting period-T orbits after perturbation based on the following elementary properties of fold maps. By smooth curve we mean a smooth arc or 1-manifold. Lemma 2.1. Let f : (R 2 , q) → (R 2 , p) be a fold singularity with singular set S(f ) (a smooth curve through q) in R 2 ; then U (f ) := f S(f ) is a smooth curve (fold image curve) through p in R 2 . Let λ be a smooth curve in R 2 passing through q, and suppose the tangent to λ at q does not coincide with the kernel of the derivative Df (q) : R 2 → R 2 . Then f (λ) is a smooth curve through p and is tangent to U (f ).

These facts can easily be deduced from the normal form (x, y) → (x, y 2 ) for a fold map (R 2 , 0) → (R 2 , 0).

Single Impact Period T Orbits

Suppose now that p ∈ Π represents a single-impact T -periodic orbit, so that p is a fixed point of the dynamical map

G c : Π → Π, that is p = G c (p) = R • φ c • F c (p).
(2.21)

This can be re-expressed as

π c (q) = R • φ c (q) (2.22)
where q = F c (p) ∈ V c . Thus, rather than seek a fixed point p of G c directly, we instead look for a point q on V c that is taken to the same point under the two maps π c and R • φ c ; the advantage here is that both of these maps are smooth.

Taking local coordinates (s, v) on V c at p 1 (here t = T + s) we find that the maps π c and φ c : V c → Π can be expressed as

π c : (s, v) → (a -1 v, v) + O(2)
(2.23)

φ c : (s, v) → (s + a -1 v, as + v) + O(2) (2.24)
for constant a = ẍ(p 1 ) = ẍ(p 0 ) > 0 . The result (2.23) is easily obtained by noting that the fold image curve is orthogonal to ( ∂xc ∂τ , ∂xc ∂v ) evaluated at p 1 and using the identity (2.32) below, while (2.24) follows from expressing V c locally near p 1 as a graph τ = τ (t, v), checking that ( ∂τ ∂t , ∂τ ∂v ) = (0, a -1 ) at p 1 and again using (2.32). Consequently, in order for (2.22) to be satisfied by a point q ∈ V c near p 1 the coordinates (s, v) of q must satisfy

a -1 v = s + a -1 v + O(2) (2.25) v = -ras -rv + O(2) (2.26)
which for small (s, v) holds only at the origin. Thus we already have the following simple but otherwise not immediately obvious result:

Proposition 2.5. Single-impact period-T orbits cannot exist arbitrarily close to a period-T grazing orbit with nonzero acceleration at the graze.

A geometric view of this fact is as follows. Equation (2.26) represents (near p 1 ) a smooth curve Γ through p 1 , namely the locus of those points q ∈ V c such that the v-coordinates of π(q) and R • φ(q) coincide. Let Γ = π c (Γ ) and

Γ = R • φ c (Γ ) c ′ V ′′ R v s v ϕ + T σ v τ = τ + σ 0 Π Π σ Γ H c c π Γ Γ ′ mod T

Fig. 2.2

The images Γ and Γ of Γ under the two fold maps πc and R • φc modulo T be the respective images of Γ : then π c (q) = R • φ c (q) precisely when the curves Γ and Γ intersect at that point. Now since Γ is tangent neither to the kernel of Dπ c (p 1 ) (the s-axis) nor the kernel of Dφ c at p 1 (which is spanned by (a, -1)) it follows from Lemma 2.1 that Γ and Γ are tangent to the fold image curves of π c and R • φ c , respectively. However, (2.23) and (2.24) show that these fold curves are tangent to the directions (1, a) and (1, -ra), respectively, and so with τ reduced modulo T these curves clearly intersect only at (τ, v) = (τ 0 , 0). See Fig. 2

.2.

From this picture it is now easy to visualize the effect of a smooth perturbation to the system, parametrized by a scalar parameter μ. We suppose that f and g in (2.1) are replaced by f μ and g μ (with c also possibly replaced by c μ ) taking their original values at μ = 0. The stability of fold singularities ensures that for sufficiently small |μ| there exists a unique point p 1 (μ) on V c (μ), varying smoothly with μ, which is a fold singularity point for both π c , φ c : V c (μ) → Π with c = c μ . Let k = k(μ) and l = l(μ) denote the τ -coordinate of the image of p 1 (μ) under π c , φ c respectively, so that the equations of the tangent lines to the fold image curves of π c , φ c through (k, 0), (l, 0) are

v = m(τ -k) (2.27) v = n(τ -l) (2.28)
with slopes m = m(μ) and n = n(μ) that are close to a and -ra, respectively, for small |μ|. It is elementary to check that the v coordinate of the point of intersection of the two lines (2.27), (2.28) is positive (so that a single impact T -periodic orbit occurs locally) if and only if kl has the same sign as

( 1 n -1 m ). Since for μ = 0 we have ( 1 n -1 m ) = -a -1 (r -1 +1
) < 0, the criterion becomes simply that l-k > 0 for sufficiently small |μ|. In the present context by construction k = τ and l = τ +t-T and so we arrive at the following result. Proposition 2.6. A single-impact T -periodic orbit occurs for sufficiently small |μ| = 0 if and only if the t-coordinate t 1 (μ) of p 1 (μ) is greater than T .

In other words, there is a T -periodic impacting orbit provided that the unique nearby initially grazing orbit that grazes a second time takes longer than time T to do so.

From the definition of p 1 (μ) as the unique point close to p 1 which satisfies x c = ẋc = 0 as well as v = 0 we find from implicit differentiation that

t 1 (0) = -a -1 ∂xc ∂τ -1 det W (2.29)
where W denotes the matrix ∂xc ∂τ ∂xc ∂μ ∂ ẋc ∂τ ∂ ẋc ∂μ evaluated at p 1 (0). If we assume t 1 (0) = 0 then since a > 0 the necessary and sufficient condition for the existence of a single-impact T -periodic orbit close to grazing becomes explicitly:

( ∂xc ∂τ ) -1 det W μ < 0. (2.30)
At first sight this appears distinct from Nordmark's criterion in Proposition 2.1 which involves derivatives with respect to v (the terms γ, δ) rather than with respect to τ . The two are reconciled, however, using the fact that, regarding the system (without the obstacle) as generating a flow {Ψ t } in (x, ẋ, τ )-space, the Jacobian matrix DΨ t (c, v, τ ) takes the initial tangent vector (v, a, 1) t to the time-t tangent vector ( ẋ, ẍ, 1) t . For a grazing periodic orbit of period T this means that (v, a, 1) t is an eigenvector of DΨ T (c, 0, τ) with eigenvalue 1, giving the identities ∂xc ∂v a + ∂xc ∂τ = 0 (2.31)

∂ ẋc ∂v a + ∂ ẋc ∂τ = a (2.32)
with derivatives evaluated at p 0 (or p 1 ). Using these we find that

det W = -a((1 -δ)κ + γλ) = -aψ
in the notation of Sect. 2.1.1. Thus the criterion (2.30) becomes γψ < 0 just as in Proposition 2.1.

Observe that here we have not required the condition det(I -A) = 0, and so (in contrast to [START_REF] Nordmark | Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators[END_REF] where this is heavily used) we have no information about the persistence of a non-impacting T -periodic orbit for |μ| = 0.

Single Impact 2T -Periodic Orbits

We next apply the impact surface geometry to the criterion for persistence of a 2T -periodic single-impact orbit. Clearly there are two criteria that together are necessary and sufficient for persistence of a 2T -periodic single-impact orbit, namely:

1. The criterion analogous to (2.30) at p 2 (0) rather than p 1 (0) 2. The condition that V c does not intervene between p 0 (μ) and p 2 (μ).

To formulate the second criterion, we need information not only about the positions of fold image curves in Π but also about which side of such a curve is the one that lies in the image of the fold map. From the normal form (x, y) → (x, y 2 ) for a fold map we deduce that a coordinate-independent description is as follows.

Lemma 2.2. Let f : (R 2 , q) → (R 2 , p) be a fold singularity as in Lemma 2.1. Let u ∈ K = ker DF (q) be a nonzero vector. The restriction f |(q + K) has the form

f : q + ξu → p + ξ 2 w + O(ξ 3 )
for a nonzero vector w ∈ R 2 that does not lie in the range of Df (q). Thus p + ηw lies in the range of f for all sufficiently small η > 0.

In our case we have that at every point p n (n > 0) the kernel K π of Dπ c (p n ) is the s-axis while the kernel K φ of Dφ c (p n ) is spanned by the vector (1, -a), using local coordinates (s, v) on V c near p n as in (2.24). Moreover, representing V c locally as the graph of a smooth function τ = τ (s, v) we find by implicit differentiation of (2.20) that at each such p n ∂ 2 τ ∂t 2 = -a ∂xc ∂τ -1 = ∂xc ∂v -1 =: γ -1 n using (2.31), and also (with a little more effort) that that is (since a number has the same sign as its inverse)

(a, -1) • D 2 φ c |K φ : ξ(-1, a) → -ξ 2 aγ -1 n det A n = -kγ -1 n for k > 0,
η((1δ)κ + γλ)μ < 0.

(2.37)

For the criterion (1) the counterpart of (2.30) at p 2 is ∂xc ∂τ

-1 2 det W 2 μ < 0,
where W 2 is the matrix (2.1.3) evaluated at p 2 . Using the easily verified fact that W 2 = (I + A)W , we may rewrite this as

η -1 γ -1 Δ + ((1 -δ)κ + γλ)μ < 0.
From (2.37) we see that the additional condition for criterion (1) to hold is that γΔ + > 0, and so we recover Proposition 2.2.

Conclusion

We have shown how the impact surface approach to the geometry of grazing bifurcation allows an easy derivation of criteria for the existence of single-impact T -periodic orbits or 2T -periodic orbits close to grazing for a T -periodically forced impact oscillator. As in the study by Nordmark [13], the methods can naturally be extended to nT -periodic orbits with m impacts. However, a further advantage of the impact surface description is that the tools of singularity theory can be applied in this context to study multiparameter bifurcation of periodic orbits close to degenerate grazing (a = 0) where the geometry of the relevant singularities (described in [START_REF] Chillingworth | Discontinuity geometry for an impact oscillator[END_REF]) is more complicated than the fold geometry that arises here. We aim to pursue this study in further work.

Fig. 2 . 1

 21 Fig. 2.1 Two fold maps πc, φc : V c → Π. The point p 1 ∈ V c is a singular (fold) point for both maps

  where the dot denotes scalar product. Using this information, sketches of the images of the fold maps π c , φ c close to p 2 and the map π c close to p 1 quickly indicate that the condition for V c not to intervene between p 0 (μ) and p 2 (μ) is close to p n = p n (0), from which we find with the obvious notation and assuming as before that γ and η (hence also γ 2 by (2.15)) are nonzero. The criterion (2.33) therefore becomes

	t 1 (μ) -t 2 (μ) γ > 0 ∂τ 0 where γ = γ 1 by definition. ∂ ẋc ∂τ a τ n (0) t n (0) = -∂xc ∂μ ∂ ẋc ∂μ with partial derivatives evaluated at p n , giving τ n (0) = -∂xc ∂τ -1 n κ n = a -1 γ -1 n κ n κμ > γ -1 2 γκ 2 μ. From the identity κ 2 = (α + 1)κ + γλ this becomes The point p ∂xc κμ > η -1 ((α + 1)κ + γλ)μ,	(2.33) (2.34) (2.35) (2.36)

n (μ) = (τ n (μ), 0; t n (μ)) is obtained by solving x c (τ, 0; t, μ) = ẋc (τ, 0; t, μ) = 0

The vector B differs from that in[START_REF] Nordmark | Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators[END_REF] by a positive scalar multiple.