Elías Cueto
email: ecueto@unizar.es

Francisco Chinesta

Meshless methods for the simulation of material forming A review

Keywords: Numerical simulation, meshless methods, forming processes

Meshless methods, that appeared in the early nineties, constitute nowadays an appealing method for the simulation of forming processes. In this review we revisit the basic ingredients of the most common of such methods, by analyzing their theoretical foundations, applicability and limitations, and give some examples of performance to show the wide variety of situations in which they can be employed.

Introduction

Although there are some examples of meshless methods dating back to the late seventies [57], the strong development of meshless methods came after the little revolution provoked by the seminal paper of Villon and coworkers on the so-called "diffuse element method" (DEM) [START_REF] Nayroles | Generalizing the finite element method: Diffuse approximation and diffuse elements[END_REF] and the popularity given to them by some modifications introduced on it by Belytschko and coworkers to create the "Element Free Galerkin" (EFG) method [START_REF] Belytschko | Crack propagation by element free galerkin methods[END_REF][START_REF] Belytschko | Element-Free Galerkin Methods[END_REF] and by W. K. Liu to create the somehow equivalent Reproducing Kernel Particle Method (RKPM) [START_REF] Liu | Wavelet and multiple scale reproducing kernel methods[END_REF].

These pioneering works opened the possibility to develop numerical methods without the need for time-consuming meshing procedures, since the connectivity of the "elements" was created by the method itself in a process transparent to the user. What is even more important, meshless methods do not suffer from lack of accuracy due to mesh (or cloud) distortion, as finite elements do. This is why, some years later, a great interest was paid on them in the forming processes community [START_REF] Chinesta | Advances in Material Forming[END_REF].

After those initial years of exploration, meshless methods came into an age of maturity once their theoretical foundations were established by Oden and coworkers [START_REF] Duarte | An h-p Adaptive Method using Clouds[END_REF], and Babuška and coworkers [START_REF] Babuška | The partition of unity finite element method: Basic theory and applications[END_REF][START_REF] Babuška | The partition of unity method[END_REF]. The flexibility, nice properties such as arbitrary degree of consistency or continuity, the neverbefore obtained spectacular simulations obtained with these methods, and also, why not, the understanding they provide on the sound properties of finite elements, provoked a decade or more of strong popularity of meshless methods and very active research interest on the numerical community. A plethora of modifications or (very often slight) improvements over existing methods gave rise to an almost endless list of different names for methods that share some common characteristics (notably, the lack of sensitivity to mesh distortion) but that nowadays remains a difficulty for those coming for the first time to the field.

After these years of maturity, only some methods came into play with really different characteristics to the ones mentioned before. Among them we can cite the Natural Element Method (NEM) [101], or the Maximum Entropy (MaxEnt) methods [START_REF] Sukumar | Construction of polygonal interpolants: a maximum entropy approach[END_REF][START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF]. They respond to some of the most important criticisms of meshless methods, namely the lack of true interpolation along the boundary of the domain, which leads to difficulties on the imposition of essential boundary conditions, and (although only partially in the case of MaxEnt) numerical integration errors.

In this paper a review is made on the use of meshless methods, with a particular emphasis on their usage in the framework of material forming, irrespective of the particular process considered. There have been some prior reviews on the field of meshless methods and, although some of them are nowadays somewhat old, we recommend the interested reader to consult [START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF][START_REF] Belytschko | Encyclopaedia of Computational Mechanics[END_REF][START_REF] Li | Meshfree and particle methods and their applications[END_REF][START_REF] Nguyen | Meshless methods: A review and computer implementation aspects[END_REF], to name a few. There are also some books available on the topic, such as [START_REF] Liu | Mesh Free Methods: Moving Beyond the finite element method[END_REF], some of them specifically devoted to aspects related to forming processes [START_REF] Chinesta | La méthode des éléments naturels dans la simulation des structures et des procédés[END_REF][START_REF] Chinesta | Natural element method for the simulation of structures and processes[END_REF].

During these years, meshless methods have successfully been applied to the simulation of a variety of forming processes, involving both solid and fluids, but with the main characteristic of employing, in general, an updated Lagrangian perspective for the description of the equations of motion. This is perhaps the most relevant novelty that meshless methods introduced, when viewed form the forming process point of view: they successfully overcome the traditional difficulty related to mesh distortion (or numerical diffusion, if prefer to employ extensive remeshings, as in many commercial codes in the field).

In this paper a review is made of the most relevant meshless (or meshfree) approaches to the field of material forming. First, a brief overview of the theoretical aspects is made in Section 2. In it, common aspects of all meshless methods are reviewed, as well as those characteristics that truly differentiate them, far beyond the long list of different names for methods that differ only in very subtle details. These include the different approximation schemes, the numerical integration schemes available to perform quadrature of the weak form of the equations, and the imposition of essential boundary conditions. These are some of the most prominent aspects of meshless methods in today's literature.

To provide the reader with a clear picture on where do we stand on the field, Section 3 includes some interesting examples of application of meshless methods that explore the ability of these methods to be applied under general conditions for the simulation of material forming processes.

Theoretical foundations of meshless methods

In this section an overview of the theoretical foundations of meshless methods is made. Particular attention is paid to the three most relevant aspects of any of such methods, namely, its approximation scheme (strictly related to its degree of continuity and consistency), its numerical integration scheme 1 and aspects related to the imposition of essential boundary conditions. These aspects are here judged as the most relevant ones for the field of forming processes. For instance, the degree of consistency of each method (the order of the polynomial they are able to reproduce exactly) is of utmost importance in the development of stable (LBB-compliant) approximations when dealing with incompressible media, as in plasticity or fluid mechanics. The degree of continuity is generally not so important, but having smooth (differentiable) approximations is essential when developing shell formulations for sheet metal forming, for instance. In turn, numerical integration is a well known source of error when dealing with meshless approximations, due to their inherent nonpolynomial character, and this is so irrespective of the particular application. Finally, many physical phenomena occur near the boundaries of the domain (contact, friction, merging flows, among others), so reproducing accurately the essential field of the problem in these regions is again a key aspect. That is why the imposition of essential boundary conditions or, more generally, the accurate interpolation of the displacement (velocity) field along the boundary is often required, and not always achieved, for meshless methods.

Meshless approximations

Probably the best form to integrate a meshless method within an existing finite element solver is to think of a finite element as the set of a particular approximation scheme (in this case, an interpolation scheme formed by piecewise 1 We focus here on methods based on the weak form of the problem, although many methods exist that are based upon the strong form, and that utilize a collocation approach, see [START_REF] Belytschko | Encyclopaedia of Computational Mechanics[END_REF]. However, their applicability to the simulation of forming processes is somewhat lower, due to aspects such as the description occurring in the vicinity of the boundary of the domain -contact, friction-and therefore are not included here. We refer the interested reader to the before mentioned reviews or books on the field for detailed explanations. polynomials) and an integration cell where Gauss quadrature is performed (again, in the case of finite elements, the intersection of the supports of its nodal shape functions). By combining different approximation schemes and different quadrature schemes, as will be seen, different meshless methods can be obtained. We review (not exhaustively) in this section the most important ones in a historic sense. In the review, moving least squares, natural neighbor, and local maximum entropy approximations are reviewed. Although not the only ones in the literature, these three approximation schemes somewhat represent three different ages in the development of meshless methods, for reasons that will be clear soon.

Moving Least Squares approximation

Although not the first in a strict historical sense, the responsible for the popularity of meshless methods is the so-called Diffuse Element Method (DEM) [START_REF] Nayroles | Generalizing the finite element method: Diffuse approximation and diffuse elements[END_REF] and, notably, the Element-Free Galerkin method (EFG) [START_REF] Belytschko | Crack propagation by element free galerkin methods[END_REF][START_REF] Belytschko | Element-Free Galerkin Methods[END_REF]. Both are based upon the approximation of the essential field of the problem (here, assume the displacement) through a moving least squares approximation. We do not include here the Smooth Particle Hydrodynamics method (SPH) [START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF]57,[START_REF] Li | Meshfree and particle methods and their applications[END_REF], since it lacks of some very basic properties such as linear consistency, nor the Reproducing Kernel Particle Method (RKPM) [START_REF] Liu | Wavelet and multiple scale reproducing kernel methods[END_REF][START_REF] Liu | Reproducing kernel particle methods[END_REF][START_REF] Gosz | Admissible approximations for essential boundary conditions in the reproducing kernel particle method[END_REF], that has been shown to be equivalent, despite its very different origin, to the EFG [START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF].

In general, not only in EFGM, the domain Ω is discretized by cloud of nodes, rather than a mesh, as seen in Fig. 1. In EFG methods, the essential field (usually, the displacement or velocity fields) is approximated as

u h (x) = m i=1 p i (x)a i (x) ≡ p T (x)a(x), (1)
where m refers to number of terms in the basis, p i (x) constitutes a polynomial basis up the desired order and, finally, a i (x) are the coefficients to determine, that notably depend on x.

As a polynomial basis, in one dimension, the linear one is often used

p T = (1, x),
or a quadratic one p T = (1, x, x 2 ), which in two dimensions reads

p T = (1, x, y, x 2 , xy, y 2 ).
The approximation given by Eq. ( 1) can be made local through 1 Covering of the domain Ω by the shape functions associated to each node, whose support is denoted by Ω I .

u h (x, x) = m i=1 p i (x)a i (x) ≡ p T (x)a(x), WI G WI G Fig.
where coefficients a i are obtained by minimization of a functional composed by the difference between the local approximation to the sought function and the essential field itself, in a least squares sense, i.e., by minimizing the following quadratic functional

J = I w(x-x I )(u h (x, x I )-u(x I )) 2 = I w(x-x I ) i p i (x I )a i (x) -u I 2 ,
where I = 1, . . . , n represents the number of nodes in the model.

Weighting, and the corresponding local character of the approximation, is given by the function w(x -x I ).This function is often a gaussian or a cubic spline. We can then re-formulate the functional as

J = (P a -u) T W (x)(P a -u),
where, as usual,

u T = (u 1 , u 2 , ..., u n ),
(in a one-dimensional case) and

P P =      p 1 (x 1 ) p 2 (x 1 ) • • • p m (x 1 ) p 1 (x 2 ) p 2 (x 2 ) • • • p m (x 2 )
. . . . . . . . . . . .

p 1 (x n ) p 2 (x n ) • • • p m (x n )      , where, finally, W W =      w(x -x 1 ) 0 • • • 0 0 w(x -x 2 ) • • • 0 . . . . . . . . . . . . 0 0 • • • w(x -x n )      .
To determine the coefficient's value, it is then necessary to minimize

J ∂J ∂a = A(x)a(x) -B(x)u = 0.
Matrix A is called matrix of moments and has the following expression

A = P T W (x)P B = P T W (x), such that a(x) = A -1 (x)B(x)u, (2)
or, equivalently,

u h (x) = n i=1 φ k i (x)u i ,
where

φ k = [φ k 1 , φ k 2 , ..., φ k n ] = p T (x)A -1 (x)B(x)
, and where k denotes the order of the approximation. The resulting shape function, for the quadratic case, is depicted in Fig. 2 for the interval [0, 1] discretized with eleven nodes. Different support sizes for W are considered, from r = 2h, with h the nodal spacing, to r = 4h.

Duarte and Oden [START_REF] Duarte | An h-p Adaptive Method using Clouds[END_REF] studied the EFG method and established for the first time one key characteristic of the method: that the shape functions constitute a partition of unity and that, despite the order of the approximation (there called intrinsic), there is the possibility of establishing and extrinsic enrichment of these functions so as to be able to make them reproduce a polynomial or arbitrary order. Melenk and Babuška generalized this approach by defining the so-called Partition of Unity method [START_REF] Babuška | The partition of unity finite element method: Basic theory and applications[END_REF][START_REF] Babuška | The partition of unity method[END_REF]: where β ji represent the new unknowns of the problem (additional degrees of freedom per node) and p j represent a basis including monomials up to a certain degree (here, care must be taken so as not to produce linear dependencies with the basis φ i , see details in [START_REF] Babuška | Meshless and Generalized Finite Element Methods: A Survey of Some Major Results[END_REF]).

u h (x) = i φ i (x) j β ji p j (x
In general, EFG methods acquired a great popularity for some years, but still present some notable drawbacks. Some of them are common for most meshless methods, others are particular of EFG. Among them we can cite the imposition of essential boundary conditions (due to the influence of interior nodes on the boundary, see Fig. 1), and the errors related to numerical integration (due to the use of quadrature cells not conforming to shape functions' support and the use of non-polynomial shape functions, see Eq. ( 2)), that will be discussed specifically in Sections 2.2 and 2.3, respectively.

Natural neighbor approximation

In the quest for a method free of errors in the interpolation of the essential variable along the boundary, the Natural Element Method (NEM) [101, 100, 102] was the first successful attempt. The NEM was originally a Galerkin method in which interpolation was achieved through natural neighbor (NN) methods [START_REF] Sibson | A Vector Identity for the Dirichlet Tesselation[END_REF][START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF][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][START_REF] Belikov | Non-Sibsonian interpolation on arbitrary system of points in Euclidean space and adaptive generating isolines algorithm[END_REF].

NEM, as was the case in EFGM, also construct the connectivity of each integration point in a process transparent to the user (and thus the appearance of no need of any mesh), but relies on the concept of Delaunay triangulations to do it. Although Delaunay triangulations are widely used in FE technology to construct meshes (it can be demonstrated that it is the best possible mesh in two dimensions), here the advantage is not the lack of any mesh, but the good accuracy provided by the NEM despite the quality or distortion of this triangulation, as proved in [101] for the first time and also in [START_REF] Cueto | Overview and recent advances in Natural Neighbour Galerkin methods[END_REF][START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF][START_REF] Cueto | Modeling three-dimensional piece-wise homogeneous domains using the α-shape based Natural Element Method[END_REF][START_REF] Alfaro | A study on the performance of natural neighbour-based galerkin methods[END_REF], among other references.

The Delaunay triangulation

D [40] of a cloud of nodes X = {x 1 , x 2 , . . . , x N } ⊂ R d , d = 2, 3,
is the unique triangulation of the cloud that satisfies the so-called circumcircle criterion, i.e., no node of the cloud lies within the circumcircle of any triangle (see Fig. 3).

The dual structure of the Delaunay triangulation is the Voronoi diagram. It is composed by a tessellation of the space into cells of the form:

T I = {x ∈ R n : d(x, x I ) < d(x, x J ) ∀ J = I},
where T I represents the Voronoi cell and d(•, •) the Euclidean distance. The simplest interpolation scheme that can be constructed on top of this geometrical construction is the so-called Thiessen interpolation, a piecewise constant interpolation within each Voronoi cell [104], which is therefore of continuity C -1 and that has been used for mixed velocity-pressure approximations in [START_REF] Gonzalez | Volumetric locking in Natural Neighbour Galerkin methods[END_REF][START_REF] Gonzalez | A natural element updated lagrangian strategy for free-surface fluid dynamics[END_REF], among others.

But the undoubtedly most popular natural neighbor interplant is due to Sibson [START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF]. If we define the second-order Voronoi cell as

T IJ = {x ∈ R n : d(x, x I ) < d(x, x J ) < d(x, x K ) ∀ K = J; ∀ K = I},
then, the Sibsonian shape function is defined as (see Fig. 4):

φ sib I (x) = κ I (x) κ(x) ,
where κ(x) and κ I (x) represents the Lebesgue measure of the cells T x and T xI , respectively. Thus defined, natural neighbor interpolation has some remarkable properties, if compared to EFG shape functions. For instance, NN shape functions are smooth (C 1 ) everywhere except from the nodes, where they are simply continuous. A simpler form to calculate NN interpolation (termed Laplace interpolation) that involves area computations instead of volumes (but produces less smooth shape functions) has been proposed in [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][START_REF] Belikov | Non-Sibsonian interpolation on arbitrary system of points in Euclidean space and adaptive generating isolines algorithm[END_REF]. Also, an alternative definition has also been proposed by Hiyoshi and Sugihara to achieve any degree of continuity in [START_REF] Hiyoshi | Two generalizations of an interpolant based on Voronoi diagrams[END_REF]. See Fig. 5 for a comparison of different shape functions on a regular grid of nodes. All the before mentioned interpolation schemes posses linear consistency, but when needed (as in the case of plasticity, where incompressibility restraints needs for the use of LBB-compliant mixed interpolations) higher-order NN interpolations can also be achieved [START_REF] Gonzalez | Higher-order natural element methods: Towards an isogeometric meshless method[END_REF]. Even a combination of NN with bubble shape functions seems to give good result for incompressible media [108].

In addition to these interesting properties, in the pioneer work of Sukumar it was reported they Sibson interpolation was interplant along convex boundaries (which is in sharp contrast to EFGM, for instance). Albeit in non-convex domains errors of about 2% were reported due to the lack of true interpolation. Later on, in [102], Sukumar claimed the interpolatory character of Laplace functions along the boundary, although in [START_REF] Cueto | On the imposition of essential boundary conditions in Natural Neighbour Galerkin methods[END_REF] some counter-examples were found that demonstrate that in concave domains this interpolation may not be true and proposed some nodal spacing restrictions so as to ensure a proper imposition of essential boundary conditions. For a deeper analysis of the issue of imposition of essential boundary conditions we refer the reader to Section 2.2.

Despite these interesting properties of NN interpolation (and others explained in Section 2.2), the main drawback of NEM is perhaps its high computational cost, especially for Sibson interpolation. In [START_REF] Alfaro | A study on the performance of natural neighbour-based galerkin methods[END_REF] a deep analysis of the computational cost of several meshless methods was accomplished. In it, it was shown that mesh distortion could lead to important inaccuracies when dealing with finite elements in practical applications, while Sibson interpolation is several orders of magnitude heavier to compute than traditional piecewise polynomial shape functions for finite elements. However, in non-linear computations, while frequent Newton-Raphson iterations are needed, the relative cost of shape function computation is obscured by the cost of updating tangent stiffness matrices.

Y Z X (a) Y Z X (b) X Y Z (c) X Y Z (d)

Local maximum entropy approximation

In the quest for the "perfect" meshless method, local maximum entropy (hereafter max-ent) methods are maybe the last to come into play. Sukumar seems to have been the first in employing global max-ent methods to solve PDEs [START_REF] Sukumar | Construction of polygonal interpolants: a maximum entropy approach[END_REF]. Max-ent approximations provide smoothness (not easily attainable by NEM), interpolation on the boundary (not easy to obtain with EFG), but are non-local in nature. That is why the original work of Sukumar proposed to use this type of interpolation on polygonally-shaped cells, thus avoiding full matrices.

More recently, Arroyo and Ortiz [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] proposed a new form of local max-ent approximation, by viewing it as a problem of statistical inference in which the in ear consistency requirement is set as a restriction to the problem. A new parameter controls the support of the shape function, thus making unnecessary to construct it over polygons.

To see the form of this new way of constructing meshless methods, consider as usual a set of nodes X = {x 1 , x 2 , . . . , x N } ⊂ R d . Let u : convX → R be a function whose values {u I ; I = 1, . . . , N } are known on the node set. "conv" represents here the convex hull of the node set. Consider an approximation of the form

u h (x) = N I=1 φ I (x)u I ,
where the functions φ I : convX → R represent the shape or basis functions. These functions are required, to be useful in the solution of second-order PDEs, to satisfy the zeroth and first-order consistency conditions:

N I=1 φ I (x) = 1, ∀x ∈ convX, ( 3a
) N I=1 φ I (x)x I = x, ∀x ∈ convX. (3b)
If these shape functions are, in addition, non-negative (φ I (x) ≥ 0 ∀x ∈ convX), then, the approximation scheme given by Eq. (2.1.3) is referred to as a convex combination, see for instance [START_REF] Farin | Curves and surfaces for CAGD[END_REF]. Since we look for positive functions φ I , the max-ent rationale is based upon the consideration of such functions as probability measures [START_REF] Sukumar | Construction of polygonal interpolants: a maximum entropy approach[END_REF]. In this sense, the Shannon entropy of a discrete probability distribution is given by:

H(φ) = - N I=1 φ I ln φ I .
Proceeding in this way, the basis function value φ I (x) is viewed as the probability of influence of a node I at a position x [START_REF] Sukumar | Construction of polygonal interpolants: a maximum entropy approach[END_REF]. The problem of approximating a function from scattered data can thus be viewed as a problem of statistical inference. Following [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF], the optimal, or least biased, convex approximation scheme (at least from the information-theoretical point of view) is the solution

of the problem maximize H(φ) = - N I=1 φ I ln φ I , (4)
subject to φ I ≥ 0, I = 1, . . . , N,

N I=1 φ I = 1, N I=1 φ I x I = x.
Proofs of the existence and uniqueness of the solution to this problem are given in [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF].

While approximations obtained after solving the problem 4 are global by definition (and thus would produce full matrices when applied in a Galerkin framework), Arroyo and Ortiz proposed a local max-ent approximation by adding spatial correlation to the problem given by Eq. ( 4). In this way, the width of the shape function φ I can be defined [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] as

w(φ I ) = Ω φ I (x)|x -x I | 2 dx,
which is equivalent to the second moment of φ I about x I . The most local approximation is that which minimizes

W (φ) = N I=1 w(φ I ) = Ω N I=1 φ I (x)|x -x I | 2 dx,
subject to the constraints given by Eqs. (3a), (3b) and the positivity restraint.

The new problem

For fixed x minimize U (x, φ)

= N I=1 φ I |x -x I | 2 (5)
subject to φ I ≥ 0, I = 1, . . . , N

N I=1 φ I = 1, N I=1 φ I x I = x
has solutions if and only if x belongs to the convex hull of the set of points [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF]. If these points are in general position, then the problem (5) has unique solution, corresponding to the piecewise affine shape functions supported by the unique Delaunay triangulation associated with the node set X (see [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] and references therein for the proof of this assertion).

Arroyo and Ortiz [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] found an elegant solution to the problem of finding a local approximation satisfying all the interesting properties of a (global) Maximum Entropy approximation by seeking a compromise between problems (4) and ( 5):

For fixed x minimize f β (x, φ) ≡ βU (x, φ) -H(φ) (6) subject to φ I ≥ 0, I = 1, . . . , N, N I=1 φ I = 1, N I=1 φ I x I = x.
Proofs of the existence and uniqueness of problem ( 6) are also given in the before mentioned reference. Note that the evaluation of the approximation (6) does not require the solution of this problem. It is enough to solve an unconstrained minimization problem that arises from the dual form of the problem [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF]. The calculation of the shape function derivatives is also explicit, see [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF].

In general, max-ent shape functions provide with a smooth approximation of the essential field, but the max-ent shape functions recover the piece-wise linear polynomials over the Delaunay triangulation of the point set if the parameter β tends to infinity (see Fig. 6). Thus defined, the method possesses linear consistency, very much like natural neighbor interpolation. It is possible, however, to obtain a (very complex) quadratic form [START_REF] Cyron | Smooth, second order, non-negative meshfree approximants selected by maximum entropy[END_REF] or, by extending the algorithm initially designed for NEM in [START_REF] Gonzalez | Higher-order natural element methods: Towards an isogeometric meshless method[END_REF], to obtain any degree of consistency also in a max-ent approach [START_REF] Gonzlez | A higher order method based on local maximum entropy approximation[END_REF]. The computational cost of max-ent approaches is considerably less than that of NEM (and higher than the corresponding FEM approach), making it an appealing candidate for the simulation of complex forming processes [START_REF] Ortiz | Maximum-entropy meshfree method for incompressible media problems[END_REF][START_REF] Gonzalez | A high order method using max-ent approximation schemes[END_REF]. We will compare some results in Section 3 to see the relative importance of this assertion in practical implementations.

Imposition of essential boundary conditions

As introduced before, the imposition of essential boundary conditions, something very natural in finite elements, has been a kind of nightmare for meshless methods. See the excellent analysis by Huerta and coworkers in [START_REF] Fernández-Méndez | Imposing essential boundary conditions in meshfree methods[END_REF]. Other works on the topic include [START_REF] Gosz | Admissible approximations for essential boundary conditions in the reproducing kernel particle method[END_REF][START_REF] Guenther | Implementation of boundary conditions for meshless methods[END_REF]. This is so since, due to the inherent unstructured connectivity between nodes in the model, interior nodes could eventually influence the result on the boundary, see Fig. 7. In general, EFG or RKP methods, whose shape function's support is normally circular or rectangular, present this problem. Also very frequently shape functions do not verify Kronecker delta properties (i.e., the approximated function does not pass through nodal values), but this is normally easier to overcome by several methods. One of the earliest methods to impose essential boundary conditions in EFG methods was simply to couple them with a strip of finite elements along the boundary [START_REF] Krongauz | Enforcement of essential boundary conditions in meshless approximations using finite elements[END_REF][START_REF] Fernandez-Mendez | Coupling finite elements and particles for adaptivity: An application to consistently stabilized convection-diffusion[END_REF][START_REF] Huerta | Enrichment and coupling of the finite element and meshless methods[END_REF]. But obviously this method somewhat eliminates the meshless character of the approximation, although it is very simple to implement. Other techniques include constrained variational principles [START_REF] Gavete | Smoothed Particle Hydrodynamics: Theory and applications to non-spherical Stars[END_REF], penalty formulations, Lagrange multipliers, and many others, but in general focused in the lack of Kronecker delta condition [START_REF] Chen | New boundary condition treatments in meshfree computation of contact problems[END_REF] and therefore did not attain a true interpolation along the boundary.

If we restrict ourselves to the case of NEM, it was assumed that the Sibson approach was interplant along convex boundaries. In [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF] it was demonstrated that through the use of he concept of α-shapes a true interpolation could be achieved even in non-convex boundaries. α-shapes are particular instances of the so-called shape constructors. Shape constructors are geometrical techniques that enable to find the shape of the domain, described solely by a cloud of nodes, at each time step. α-shapes [START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF] have been employed in a number of previous works involving free surface flows, see for instance [START_REF] Idelsohn | The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF] [71] [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF] [59] [START_REF] Birknes | A particle finite element method applied to long wave run-up[END_REF] [88] [89], among others. In essence, through the definition of a parameter α that represents the level of detail up to which the geometry is to be represented, this technique allows to extract almost with no user intervention the actual geometry of the domain, see Fig. 8.

Therefore, in addition to the importance of the true interpolation achieved by this method, α-shapes (or, in general, shape constructors) provide a very efficient means to deal with nodal clouds evolving in time, as will be analyzed in Section 3 for fluid mechanics problems. Fragmentation, coalescence, merging flows, etc., can be treated with no user intervention very conveniently in this way.

In [107] a new approach was developed for the imposition of essential boundary conditions in the context of the NEM (Constrained-NEM, or C-NEM) that is based upon the usage of a visibility criterion, a concept initially developed in the context of EFGM [START_REF] Organ | Continuous meshless approximations for nonconvex bodies by diffraction and transparency[END_REF]. It can be demonstrated that, up to some differences in the implementation aspects, it produces results entirely equivalent to those obtained by employing α-shapes, but needs for an explicit description of the boundary of the domain in the form of a triangulation or planar straight line graph [106, [START_REF] Illoul | A mixed finite element/meshless natural element method for simulating rotative electromagnetic machines[END_REF]109,[START_REF] Chinesta | La méthode des éléments naturels dans la simulation des structures et des procédés[END_REF][START_REF] Chinesta | Natural element method for the simulation of structures and processes[END_REF][START_REF] Illoul | On some aspects of the CNEM implementation in 3D in order to simulate high speed machining or shearing[END_REF].

For Laplace-type of NN interpolations, it was initially claimed [102] that they were interplant even on non-convex domains, but in [START_REF] Cueto | On the imposition of essential boundary conditions in Natural Neighbour Galerkin methods[END_REF] it was proved that this interpolation was only node-wise and that some counter-examples could be found. An alternative method was proposed in this last reference to impose exactly essential boundary conditions by employing a planar straight line graph describing the boundary.

If we restrict ourselves to the case of local max-ent approximation, in [6] a demonstration is included that proves that the method possesses a weak Kronecker delta property, i.e., for convex domains, the value of shape functions associated to interior nodes vanishes on the boundary. Thus imposing the exact value of the nodal displacement with Lagrange multipliers, for instance, will suffice to verify essential boundary conditions straightforwardly, very much like in finite elements. This is in general not true for non-convex domains, however. In [START_REF] Gonzalez | Higher-order natural element methods: Towards an isogeometric meshless method[END_REF] this approach was later generalized to arbitrary order of consistency. See Fig. 9 for an example of quadratic max-ent shape function on the boundary of the domain verifying Kronecker delta property.

Numerical integration

As mentioned before, only formulations based upon weak forms of the governing equations are being covered in this review, since collocation methods are, in principle, not so well suited for the simulation of forming processes. In this framework, numerical integration plays a fundamental role for the accuracy of the results. In general, meshless formulation do not employ polynomial shape functions, and since quadrature formulas were initially designed to integrate such functions, an error is unavoidable. In addition, to perform numerical integration an integration cell should be chosen, and these do not generally conform with the support of shape functions (as is trivially the case in finite elements, where the integration cell arises naturally as the intersection of the element's shape functions support). This provokes a second source of error.

Historically, the first attempt to overcome problems related to numerical integration errors was to try to conform integration cells to shape function support, as in [START_REF] Dolbow | Numerical Integration of the Galerkin Weak Form in Meshfree Methods[END_REF][START_REF] De | Towards an efficient meshless computational technique: the method of finite spheres[END_REF][START_REF] De | The method of finite spheres with improved numerical integration[END_REF][START_REF] Gonzalez | Numerical integration in Natural Neighbour Galerkin methods[END_REF][START_REF] Atluri | A critical assesment of the truly Meshless Local Petrov-Galerkin and Local Boundary Integral Equation methods[END_REF][START_REF] Atluri | New concepts in meshless methods[END_REF], for instance. In general, none of these methods completely overcome the mentioned deficiencies, due to the presence of the second source of error mentioned earlier: the non-polynomial character of the approximation.

The most important advancement in the numerical integration of meshless methods arose with the development of the so-called Stabilized Conforming Nodal Integration (SCNI) by J. S. Chen [START_REF] Chen | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF]. In essence, this method is based on assuming a modified strain field at each node:

εh ij (x I ) = Ω ε ij (x)Φ(x; x -x I )dΩ,
where ε represents the Cauchy strain tensor and Φ is a distribution function, that is usually taken as:

Φ(x; x -x I ) = 1 A I if x ∈ Ω I 0 otherwise
with Ω I the Voronoi cell associated to node I (other approaches in the definition of the area associated to each node are equally possible) and A I the corresponding area of this cell. With this definition, the strain smoothing leads to:

εh ij (x I ) = Ω 1 2
∂u h i ∂x j + ∂u h j x i Φ(x; x -x I )dΩ
By the divergence theorem, it can be obtained

εh ij = 1 2A I Γ I (u h i n j + u h j n i )dΓ, (7)
Γ I being the boundary of the Voronoi cell associated to node I. By introducing shape functions into Eq. ( 7), a matrix expression is obtained in the form:

εh (x I ) = J∈N N (I) BJ (x I )u I ,
where N N (I) represents the set of nodes neighboring the point x I . The approximation to the weak form of the problem leads to a stiffness matrix and a force vector, in the absence of body forces, that can be expressed as:

K IJ = N P m=1 BT I (x m )C BJ (x m )V n m , f I = N nb m=1 φ(x m )t(x m )V n-1 m
, where V n m denotes the volume in dimension n of the Voronoi cell associated to node m. N nb represents the number of nodes in the natural boundary.

This nodal quadrature scheme has rendered excellent results when applied to the EFGM [START_REF] Chen | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF]. This technique is also very well suited for its application within the NEM, since most of the geometrical entities appearing in its computation (Voronoi cells, circumcenters, etc.) are in fact previously obtained in the NE shape function computation and can be easily stored with important computational saving.

Also specially noteworthy is the fact that SCNI leads to a truly nodal implementation of the Natural Element method (or whatever method it is applied to), which means that no recovery of secondary variables must be performed, nor nodal averaging, like in the traditional version of the Finite Element method or in previous non-linear versions of the NEM [START_REF] Garcia-Aznar | Simulation of bone internal remodeling by means of the α-shape based Natural Element Method[END_REF]. In Fig. 10 an example of application of this technique is shown over the geometry of a hollow cylinder [START_REF] Gonzalez | Numerical integration in Natural Neighbour Galerkin methods[END_REF]. It can be noticed how a Voronoi cell must be computed around each node and a quadrature scheme stablished on the boundary of it, usually by performing a division on triangles.

Surprisingly or not, W. Quak [START_REF] Quak | A comparative study on the performance of meshless approximations and their integration[END_REF][START_REF] Quak | Meshless methods and forming processes[END_REF] found that by applying SCNI to standard finite elements very good results were obtained in terms of accuracy for very distorted meshes. In this reference, a bending test was simulated with five different clouds of nodes, see Fig. 11. It was found that finite elements provided very competitive results, comparable to those of max-ent approximations, and much better than EFG methods. Nowadays there is a plethora of nodal integration techniques developed for finite elements that has provided them with the best characteristics of meshless methods in terms of robustness to mesh distortion [92, [START_REF] Elmer | A stable, meshfree, nodal integration method for nearly incompressible solids[END_REF][START_REF] Castellazzi | Displacement-based finite elements with nodal integration for Reissner-Mindlin plates[END_REF][START_REF] Krysl | Locking-free continuum displacement finite elements with nodal integration[END_REF].

Examples of application

Once meshless methods acquired an important maturity, their application to forming process simulation became almost straightforward, as will be seen. There have been many applications in the field of solid (bulk) forming processes, as well as in fluid mechanics, but slightly less in sheet metal forming, for reasons that will become clearer later on. In this section a revision is made on the different applications that can be envisaged for an advantageous application of meshless methods.

Bulk forming processes

Bulk forming processes is maybe the most natural application of meshless methods to the field of forming, and one of the first where they were applied is the effective strain rate. In particular, in [START_REF] Alfaro | Three-dimensional simulation of aluminium extrusion by the α-shape based natural element method[END_REF][START_REF] Alfaro | Meshless methods with application to metal forming[END_REF] it was assumed that the aluminum yield stress varied according to a Sellars-Tegart law:

σ y (d) = S m arcsinh d 1 A e Q RT 1 m with d 1 = max{d, d 0 }, (8)
and a coupled thermo-mechanic model was implemented. The key ingredient in this model was to employ an updated Lagrangian framework to describe flow kinematics, i.e., an explicit update of the geometry is achieved through

x n+1 I = x n I + v n I ∆t,
and where the velocity field v is obtained after a consistent linearization of the variational problem. Since no loss of accuracy is produced by 'mesh' distortion, this very simple approach to many forming processes produces very appealing results for many forming processes like aluminum extrusion (already mentioned, but also in [START_REF] Alfaro | Simulation of the extrusion of hollow profiles by natural element methods[END_REF][START_REF] Filice | A preliminary comparison between finite element and meshless simulations of extrusion[END_REF]), friction stir welding [START_REF] Alfaro | Numerical simulation of friction stir welding by natural element methods[END_REF], forging [START_REF] Sukumar | Meshless Methods and Partition of Unity Finite Elements[END_REF], casting [START_REF] Alfaro | Three-dimensional simulation of aluminium extrusion by the α-shape based natural element method[END_REF], laser surface coating [START_REF] Gonzalez | Natural Neighbour Strategies for the Simulation of Laser Surface Coating Processes[END_REF], machining [START_REF] Chinesta | Thermomechanical Cutting Model Discretisation. Eulerian or Lagrangian[END_REF], and many others. Many practical aspects of forming simulation were studied in [START_REF] Alfaro | A study on the performance of natural neighbour-based galerkin methods[END_REF]. In particular, aspects related to computational cost and accuracy were deeply studied. Natural element and finite element methods were compared and, noteworthy, it was found than the computational cost of NEM was higher, particularly if Sibson shape functions are used. However, since we deal with highly nonlineal applications, this computational cost is obscured by the time spent in obtaining consistent tangent stiffness matrices in the Newton-Raphson loop.

In Fig. 12(a) geometry of an extrusion test is show (only one half of the domain is simulated by applying appropriate boundary conditions). At this location, equivalent strain rate contour plots for finite element, Sibson-NEM and Laplace-NEM, respectively. Note how finite element approximation produces an artificial stress concentration near the symmetry plane, where nothing can provoke this spurious stress level, that can only be due to mesh distortion. On the contrary, NEM (both Sibson and Laplace approximations) produced excellent results, with no apparent spurious stress concentrations.

Specially noteworthy is the ability of meshless methods to accurately reproduce the process of porthole extrusion. Due to the particular geometrical conditions of the dies necessary to obtain hollow profiles, porthole extrusions needs for a simulation strategy able to take into account the process of material separation and welding through extrusion. In [START_REF] Gagliardi | NEM-FEM comparison on porthole die extrusion of AA-6082[END_REF] a deep comparison of FE and NEM techniques was performed for the porthole die extrusion of AA-6082 and compared with experimental results. In Fig. 13 a comparison is made with different die geometries in which different results can be noticed. In the model, it was assumed that a proper welding necessitated of a prior level of pressure within the die chamber. A nodal implementation of NEM allowed to track the pressure history at nodal positions, thus enabling to know if a proper welding had been achieved. In Fig. 13 blue nodes represent lack of true welding after the porthole, whereas red nodes indicate a good process. It can be noticed that only the geometry represented in the middle figure produces good quality profiles due to its particular geometry. The rightmost figure, despite its good apparent geometry, shows that no proper welding is achieved (note the blue nodes all along the geometry of the extrudate). Under internal pressure, it is expected that this welding will open again, as shown by experimental results.

Sheet metal forming

The corps of literature of meshless methods devoted to sheet metal forming is considerably smaller than that of bulk forming. In some way it is natural, since meshless methods do not employ any connectivity dictated by the user to define a "mesh". Therefore, shells, that constitute a manifold geometry in three-dimensional space, are difficult to reproduce by a method based solely on a cloud of nodes.

Among the first works devoted to sheet metal forming one can cite [105] and [START_REF] Liu | Adaptive multiple scale meshless simulation on springback analysis in sheet metal forming[END_REF]. In both a kind of solid-shell approach was employed, i.e., several nodes were placed through the thickness direction of the shell. Very few works were devoted to true shell formulations. Maybe [START_REF] Sidibe | A meshfree simulation of the draw bending of sheet metal[END_REF] can be considered as the first one, up to our knowledge. In all these references, since no connectivity is established, and in order to avoid computations of distances between nodes in secant directions (instead of in directions tangent to the manifold), a nodal spacing, h, It was not until very recently that a sound basis has been established for the proper simulation of shells with meshless methods [START_REF] Millan | Nonlinear manifold learning for meshfree finite deformation thin-shell analysis[END_REF], in this case by employing max-ent approaches, although the method is completely general. Briefly speaking, the proposed method is based upon the assumption that a shell is actually a manifold in geometrical terms, i.e., each point of a shell Under these assumptions, only the central geometry will provide an appropriate quality in the extruded profile. The rightmost geometry, despite its apparent good geometry, will provide a profile that will open once it is subjected to internal pressure.

resembles locally an Euclidean (flat) space. Therefore, in order to perform a consistent meshless analysis of shells, it is necessary to locally learn its manifold structure, by identifying its locally flat principal directions. This is done on a node-by-node basis by employing non-linear principal component analysis (PCA) techniques [103].

Fluid forming processes

Fluid forming processes involve a great variety of situations, ranging from casting [106,[START_REF] Alfaro | Meshless methods with application to metal forming[END_REF] to mould filling [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF], from Resin Transfer Moulding [START_REF] Garcia | Meshless methods with application to liquid composite molding simulation[END_REF] to spin coating processes [START_REF] Cueto | Numerical simulation of spin coating processes with carbon nanotubes suspensions[END_REF][START_REF] Cueto | Numerical simulation of spin coating processes involving functionalised carbon nanotube suspensions[END_REF][START_REF] Cueto | Rheological modeling and forming process simulation of cnt nanocomposites[END_REF], where a plethora of different models for polymers, among other materials, could be taken into account [START_REF] Chinesta | Modeling nanocomposites: from rheology to forming processes simulation[END_REF]. As in the case of bulk forming processes, meshless methods allow for an updated Lagrangian description of fluid flows easily [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF]. This is particularly noteworthy when free surface flows are present, and can be of little help if not. However, nodal implementations can help in designing efficient algorithms when variables depending to history are to be taken into account.

To exemplify how this affects the simulation of fluid flows, let us restrict to Navier-Stokes equations, for instance:

∇ • σ + ρb = ρ dv dt = ρ ∂v ∂t + v∇ • v , (9a)
∇ • v = 0, (9b)
σ = -pI + 2µD. ( 9c
)
The weak form of the problem associated to Eqs. (9a), (9b) and (9c) is:

Ω 2µD : D * dΩ - Ω pI : D * dΩ = - Ω ρb • v * dΩ + Ω ρ dv dt • v * dΩ, (10)
and

Ω ∇ • v p * dΩ = 0,
where ":" denotes the tensor product twice contracted and b the vector of volumetric forces applied to the fluid. D * represents and admissible variation of the strain rate tensor, whereas v * represents equivalently an admissible variation of the velocity. The second term in the right-hand side of Eq. ( 10) represents the inertia effects. Time discretization of this term represents the discretization of the material derivative along the nodal trajectories, which are precisely the characteristic lines related to the advection operator. Thus, assuming known the flow kinematics at time t n-1 = (n -1)∆t, meshless methods allow to proceed easily as follows, by virtue of the updated Lagrangian framework:

Ω ρ dv dt v * dΩ = Ω ρ v n (x) -v n-1 (X) ∆t v * dΩ, (11)
where X represents the position at time t n-1 occupied by the particle located at position x at present time t n , i.e.:

x = X + v n-1 (X)∆t.

So we arrive at

Ω 2µD : D * dΩ - Ω pI : D * dΩ - Ω v • v * ∆t dΩ = = - Ω ρb • v * dΩ - Ω ρ v n-1 • v * ∆t dΩ, (12)
and

Ω ∇ • v p * dΩ = 0.
where the superindex in all the variables corresponding to the current time step has been dropped for clarity.

The most difficult term in Eq. ( 12) is the second term of the right-hand side. The numerical integration of this term depends on the quadrature scheme employed [START_REF] Gonzalez | A natural element updated lagrangian strategy for free-surface fluid dynamics[END_REF].

If traditional Gauss-based quadratures on the Delaunay triangles are employed, it will be necessary to find the position at time t n-1 of the point occupying at time t n the position of the integration point ξ k (see Fig. 14):

Ω ρ v n-1 • v * ∆t dΩ = k ρ v n-1 (Ξ k ) • v * (ξ k ) ∆t ω k , (13)
where ω k represents the weights associated to integration points ξ k , and Ξ k corresponds to the position occupied at time t n-1 by the quadrature point ξ k , see Fig. 14.

v ( ) n-1 X v ( ) n x
Fig. 14 Determination of the position of quadrature points at time step t n-1 .

On the contrary, if some type of nodal integration, as in [START_REF] Gonzalez | Numerical integration in Natural Neighbour Galerkin methods[END_REF], is employed, this procedure becomes straightforward, with the only need to store nodal velocities at time step t n-1 .

In [START_REF] Idelsohn | The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF] a Lagrangian method that employs Natural Neighbor interpolation to construct the discrete form of the problem was presented. In that case, however, an implicit three-step fractional method was employed to perform the time integration. This approach needs for a stabilization if small time increments are chosen. See [START_REF] Idelsohn | The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF] for more details. In that reference, however, the method is truly a particle method, since nodes posses volume and mass associated to them.

If free-surface flows are considered, it is again of utmost importance to employ any technique for the reconstruction of the geometry of the evolving domain. While most meshless approaches do not employ any particular technique, thus provoking results depending on the shape functions' support, shape constructors have been employed in a number of works [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF][START_REF] Gonzalez | A natural element updated lagrangian strategy for free-surface fluid dynamics[END_REF][START_REF] Idelsohn | The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF][START_REF] Idelsohn | To mesh or not to mesh: that is the question[END_REF][START_REF] Cante | Flow regime analyses during the filling stage in powder metallurgy processes: experimental study and numerical modelling[END_REF][START_REF] Hartmann | A 3D Frictionless Contact Domain Method for Large Deformation Problems[END_REF]. α-shapes [START_REF] Edelsbrunner | On the shape of a set of points in the plane[END_REF][START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF][START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF] are perhaps the most widely used shape constructors in this context. However, it is well known that this technique usually fails in the presence of holes or merging flows during the simulation, unless a very fine nodal cloud is employed. In [START_REF] Galavis | Improved boundary tracking in meshless simulations of free-surface flows[END_REF] an improved α-shape method is developed that takes into account not only geometrical features of the cloud of points, but also the history of velocity field in order to anticipate merging flows, the appearance of holes, etc, to a higher degree of accuracy. This method has rendered excellent results in the simulation of free-surface flows and also in the simulation of fluid-structure interactions [START_REF] Galavis | A natural element updated lagrangian approach for modelling fluid structure interactions[END_REF].

Particularly noteworthy is the ability of meshless methods to construct node-based approaches for fluids with complex constitutive equations [START_REF] Galavis | A natural neighbour Lagrange-Galerkin method for the simulation of Newtonian and Oldroyd-B free surface flows[END_REF]. Multiscale methods arising from kinetic theory [START_REF] Oettinger | Stochastic Processes in Polymeric Fluids[END_REF] are a clear example of this. In these models, a constitutive equation of the type

σ = -pI + τ + 2ηD,
where τ represents an extra-stress contribution coming from the micro-scale. To obtain this extra-stress contribution, kinetic theory provides also the equation governing the evolution of the probability distribution function ψ at the micro-level. This equation is known as the Fokker-Planck equation:

Dψ Dt = - ∂ ∂X • {Aψ} + 1 2 ∂ ∂X ∂ ∂X : {Cψ}, (14)
where A represents a vector describing the drift exerted by the fluid on the function ψ, and C is a symmetric, positive-definite matrix accounting for brownian effects in the model. D/Dt represent material derivative. The expression, finally, that relates the obtained configuration state with its enforced state of stress is known as the Kramers formula:

τ = g(X)ψdX = g(X)
where the brackets • denote an ensemble average over all the molecular conformational space at a physical point, and g is some function of the configuration state, depending on the particular model considered.

The stochastic approach to these problems makes use of the equivalence of the Fokker-Planck Eq. ( 14) to the following Itô's stochastic differential equation [START_REF] Oettinger | Stochastic Processes in Polymeric Fluids[END_REF]:

dX = Adt + B • dW (15)
where

C = B • B T
and W represents a Wiener process. Eq. ( 15) applies along individual molecule trajectories.

Meshless methods allow to attach molecule ensembles to the nodes, allowing for the integration of Eq. ( 15) by the method of characteristics, along the nodal paths straightforwardly [START_REF] Cueto | MESHLESS STOCHASTIC SIMULATION OF MICRO-MACRO KINETIC THEORY MODELS[END_REF]. The integration is performed by means of the so-called Euler-Maruyama scheme: where n refers to the current time step and j to the individual molecule being integrated.

X j n+1 = X j n + A(X j n , t n )∆t + B(X j n , t n ) • ∆W j n (16)
For instance, in [START_REF] Cueto | MESHLESS STOCHASTIC SIMULATION OF MICRO-MACRO KINETIC THEORY MODELS[END_REF] an study is made by means of natural elements of the swelling behavior of an entangled polymer modeled by a Doi-Edwards fluid. Four snapshots of the velocity field are shown in Fig. 15. Elastic effects are notorious after the outlet of the channel. A minor loss of symmetry in the flow is noticed due mainly to the statistical noise. Despite the very low number of realizations per node, the statistical noise remains surprisingly low.

In [START_REF] Galavis | A natural neighbour Lagrange-Galerkin method for the simulation of Newtonian and Oldroyd-B free surface flows[END_REF], for instance, it is shown how meshless methods provide an approximation to the swelling flow of non-newtonian flows that improves the accuracy of existing finite element approximations.

Conclusions

Meshless methods arose in the middle nineties as a promising alternative to finite elements where the process of generation of complex meshes constitutes a major issue or where mesh distortion provoke loss of accuracy in the results. Meshless methods provide a very flexible alternative for these cases, although all the problems they presented initially needed for a very active research activity during more than a decade. Today, some twenty years after, meshless methods have overcome most of their initial limitations and constitute nowadays an appealing alternative in many fields. Noteworthy, material forming simulation is one of such fields.

During this review it has been highlighted how meshless methods are now able to provide a very competitive alternative to finite element simulation in many fields, from bulk forming to complex fluid flows (and perhaps to a lesser extent in sheet metal forming). Their main ability is perhaps to provide a very convenient way of performing updated Lagrangian descriptions of forming processes. This type of description is helpful when dealing with bulk forming, for instance, or when free-surface flows are present.

Noteworthy, all the research activity generated by meshless methods has helped to improve also the properties of finite elements. It has been shown, for instance, how finite elements with stabilized conforming nodal integration provide very competitive results at a fraction of the computational cost of meshless methods. All this research activity is today, after some 20 years, helping us to understand why finite elements, originated some sixty years ago, are still among us, probably for many years to come.

Fig. 2

 2 Fig. 2 Example of quadratic EFG shape function on a regular one-dimensional grid for different support sizes.

Fig. 3

 3 Fig. 3 Delaunay triangulation and Voronoi diagram of a cloud of points. The rightmost figure shows a degenerate situation in which four nodes lye in the same circle.

Fig. 4

 4 Fig. 4 Definition of the Natural Neighbor coordinates of a point x.

Fig. 5

 5 Fig. 5 Sibson (a), Laplace (b), and Hiyoshi-Sugihara interpolants with C 1 (c) and C 2 (c) continuity, respectively.

Fig. 6 Fig. 7

 67 Fig.6Influence of the β parameter on the resulting shape function. Functions φ I (x) for the point located at the centre of the cloud and parameters β = 0.2, 0.6, 2.0 and 4.0, respectively. Note the different supports, but also the different heights of the functions on the scale.

Fig. 8

 8 Fig. 8 Evolution of the family of α-shapes of a cloud of points representing a wave breaking on a beach. Different shapes for different α values (S α=0 or cloud of points (a), S 0.5 (b), S 1.0 (c), S 2.0 (d), S 3.0 (e) and S∞ (f)) are depicted.

Fig. 9

 9 Fig.9Interpolating max-ent quadratic function on the boundary of the convex hull of the data sites. A discontinuity on the derivative appears at the node location in this case.

Fig. 10

 10 Fig.10Division of a hollow cylinder[START_REF] Gonzalez | Numerical integration in Natural Neighbour Galerkin methods[END_REF] into Voronoi cells. Note that concave regions need for an additional triangulation to intersect Voronoi cells, defined on the convex hull of the node cloud.

Fig. 11

 11 Fig. 11 Distortion test performed with Stabilized Conforming Nodal Integration.

  [START_REF] Chen | Analysis of metal forming process based on meshless method[END_REF][START_REF] Alfaro | Three-dimensional simulation of aluminium extrusion by the α-shape based natural element method[END_REF][START_REF] Alfaro | Meshless methods with application to metal forming[END_REF]. One of the most classical formulations when dealing with processes like aluminum extrusion, for instance, is to assume that, since plastic strains are much larger than elastic ones, one can neglect them and employ the socalled flow formulation (see[112] [110] [111][START_REF] Lof | FEM simulations of the extrusion of complex thin-walled aluminium sections[END_REF], just to cite some of the first and more recent works using this assumption) to assume that aluminum behaves like a non-Newtonian fluid, i.e., ε = f (σ).By assuming a rigid-viscoplastic constitutive law, and a von Mises plasticity criterion, one arrives at σ = 2µd -pI, with µ = σ y 3d ,

Fig. 12

 12 Fig. 12 Contour plot of the second invariant of the strain rate tensor at the location indicated in (a). (b) FEM, (c) Sibson, (d) Laplace results.

Fig. 13

 13 Fig. 13 Porthole die extrusion under different die geometries. Red nodes indicate a good quality in the welding of aluminum after the porthole, whereas blue nodes indicate defects.Under these assumptions, only the central geometry will provide an appropriate quality in the extruded profile. The rightmost geometry, despite its apparent good geometry, will provide a profile that will open once it is subjected to internal pressure.

Fig. 15

 15 Fig. 15 Evolution of the velocity field in the die swelling flow of an entangled polymer.

This work has been partially supported by the Spanish Ministry of Economy and Competitiveness, through grant number CICYT-DPI2011-27778-C02-01.