# A simple error indicator for meshfree methods based on natural neighbors

Julien Yvonnet, Gérard Coffignal, David Ryckelynck, Philippe Lorong, Francisco Chinesta

► To cite this version:

Julien Yvonnet, Gérard Coffignal, David Ryckelynck, Philippe Lorong, Francisco Chinesta. A simple error indicator for meshfree methods based on natural neighbors. Computers & Structures, 2006, 84 (21), pp.1301 - 1312. 10.1016/j.compstruc.2006.04.002. hal-01004942

HAL Id: hal-01004942

https://hal.science/hal-01004942

Submitted on 15 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
A simple error indicator for meshfree methods
based on natural neighbors

J. Yvonnet<sup>a</sup>, G. Coffignal<sup>b</sup>, D. Ryckelynck<sup>b</sup>, Ph. Lorong<sup>b</sup>, F. Chinesta<sup>b</sup>

<sup>a</sup> Université de Marne-la-Vallée, LAM, 5 Bd Descartes, 77454 Marne-la-Vallée Cedex, France

<sup>b</sup> LMSP, UMR 8106 CNRS ENSAM-ESEM, 151 boulevard de l'Hôpital, F-75013 Paris, France

The main aim of this paper is the development of a refinement procedure able to operate in the context of the constrained natural element method (C-NEM). The C-NEM was proposed by the authors in a former work and its main meshless features, that allow to describe large domain changes as well as to handle fixed or moving discontinuities, were analyzed. Sometimes, in order to improve the interpolation accuracy for describing boundary layers or an anisotropic behavior, new nodes must be added, removed or repositioned. The interpolation in the vast majority of meshless techniques is free of mesh quality requirement. Thus, introduction, elimination or repositioning of nodes is a trivial task, because no geometrical restrictions exist. In this way, nodes can be added without geometrical checks in the regions where the solution must be improved (identified by using an appropriate error indicator). For this purpose, in this paper an a posteriori error indicator will be proposed and tested in some linear elastostatic problems benchmarks involving different levels of difficulty (stress concentration, solution singularities, ...) all of them with a known exact solution. The computational implementation of this error indicator is very simple, and when it is used in tandem with an efficient refinement procedure, which makes use of the meshless features of the C-NEM, provides an accurate adaptation procedure, specially appropriate in the C-NEM framework.

*Keywords: Adaptivity; Error indicators; Meshless methods; Natural element method; Constrained natural element method*

# 1. Introduction

The objective of this work is to define an efficient error indicator and the associated refinement procedure within the framework of the meshless natural neighbor Galerkin method. The finite element method (FEM), which has been widely used in many engineering problem simulations, exhibits some limitations when the mesh become highly distorted. It is well known that finite element interpolation

fails when the elements (where the interpolation is defined) becomes too distorted [2]. In order to avoid these problems, alternative approaches, known as meshless methods, have been developed. Meshless methods describe a continuum body by a finite number of nodes. Thus, the unknown fields are interpolated everywhere from the nodal values that constitute the problem degrees of freedom. The main advantage of meshless methods is the fact that the interpolation accuracy is much less affected by the nodal distribution. However, it is obvious that in any case, an appropriate nodal density is required in order to describe high gradients (boundary layers) as well as an anisotropic behavior of
the solution. For this reason, nodal adaptation is needed to
compute numerical solutions of problems governed by par-
tial differential equations.

In the framework of the finite elements, these adaptation
procedures are known as remeshing techniques. Remeshing
is required for example when the element geometry
becomes too distorted as a consequence of large domain
changes. Moreover, sometimes, in order to improve the
interpolation accuracy for describing boundary layers or
an anisotropic behavior, new nodes must be added,
removed or repositioned, and in fact this is not an easy task
because the mesh associated to the new nodal distribution
cannot contain any too distorted element. If this is not the
case, local or global remeshing is required in order to guar-
antee the geometrical quality of the mesh elements.

On the contrary, in meshless techniques, interpolation is
free of that mesh requirement. Thus, introduction, elimina-
tion or repositioning of nodes is a trivial task, because no
geometrical restrictions exist. In this way, nodes can be
added without geometrical checks in the regions where
the solution must be improved (identified by using an
appropriate error indicator). Once that the new nodes are
placed into the domain, and for problems making use of
internal variables, these variables can be initialized at those
nodes using the standard meshless interpolation. This
appealing feature of this kind of techniques simplifies sig-
nificantly the refinement procedures. Moreover, the use
of smoother shape functions (constructed without a finite
element mesh support) allows to deal with large geometri-
cal transformations [31].

Many meshless methods have been developed, including
the meshfree finite differences [21], the smooth particle
hydrodynamics (SPH) [30], the diffuse element method
(DEM) [33], the h-p clouds method [14], the element-free
Galerkin (EFG) method [4], the reproducing kernel parti-
cle method (RKPM) [27], the radial basis functions
(RBF) and the method of finite spheres [12], among many
others. Nevertheless, a troublesome task in these tech-
niques is the imposition of essential boundary conditions,
because in general the Kronecker delta property is not sat-
isfied. Despite of several techniques have been proposed to
overcome this problem [4,7,19], an appealing alternative
technique which allows to circumvent this difficulty lies
in the use of the natural neighbor interpolation (NEM)
[34,36].

The natural neighbor Galerkin method, or natural ele-
ment method (NEM), which has been considered as a
member of the meshless methods family, satisfies the Kro-
necker delta property as well as the linear consistency.
Sukumar has shown in [36] the strict linearity of the shape
functions along any convex boundary. However, on the
non-convex ones, the NEM shape functions are no more
linear. In this way, we consider a point located on a non-
convex boundary. The interpolated value of a field at this
position depends on the two neighbor nodes located on
the boundary and on other spurious nodes: some are
located in the neighborhood of this point within the

domain and others on the non-convex boundary far from
the considered point, as described in our former works [42].

A strategy to overcome this problem and to recover the linearity of the shape functions along any kind of boundary (convex or not) was proposed by Cueto et al. in [11] by invoking the concept of alpha shapes. Moreover, this technique allows to proceed without an explicit description of the domain boundary, which is automatically extracted from the cloud of nodes. However, the treatment of problems involving cracks requires an appropriate adaptation of the nodal density, which can reach infinity as the width of the crack goes to zero.

On the contrary, if the domain boundary is defined
using a CAD description, the shape functions can be com-
puted from the constrained Voronoi diagram instead of the
Voronoi diagram used in the NEM, and no more nodal
density control is required. It has been proved in our for-
mer works [42] that the use of this technique, known as
constrained natural element method (C-NEM), allows to
recover the shape functions linearity along any kind of
boundary. Despite of the fact that the constrained Voronoi
diagram makes use of a visibility criterion, no numerical
discontinuities are introduced. Thus, cracks and fixed or
moving discontinuities can be treated in a very simple
manner [42,43].

Error estimation and the associated remeshing are well established in the finite elements framework [22,3,47]. A state of the art on basic a posteriori error estimation, including limits of usual error measures, has been proposed recently in [16]. Nevertheless, the emergence of papers on these topics in the context of meshless methods is more recent, despite of its high simplicity, as previously argued, for the insertion and deletion of nodes. In the last years, some researches have been devoted for developing suitable adaptive approaches for meshless methods.

In the context of the reproducing kernel particle method (RKPM), Belytschko et al. [5], and Liu et al. [28] have proposed approaches based on the residuals. You et al. [41] have used an original approach by utilizing the reproducing kernel as a low-pass filter and the corresponding high-pass filter is used to identify the locations of high gradient and serves as an operator for error indication. In the context of the element-free Galerkin (EFG), error estimation has been widely investigated using recovery based error indicators [47]. Owing to the smooth shape functions in meshless methods, the stress recovery error does not require any projection as was needed in the C⁰ finite element methods. Chung and Belytschko [10] introduced local and global error estimates for the element-free Galerkin (EFG) method. The essence of this error estimate is to use the difference between the values of the recovered stress and those obtained directly by the EFG solution. More in-depth study about refinement procedures in EFG have been conducted by Lee and Zhou [24,25] and in the RKPM by Lu and Chen in [29]. Nevertheless, some difficulties appear, because the quality of the effectivity index depends on the number of nodes contained in the domain of
influence used in the computation and in the projection stages [10]. Moreover, the difficulties associated to the adaptation of the domain of influence during refinement are discussed in [24].

The main aim of the present work is the development of a refinement procedure able to operate in the context of the C-NEM, which includes an a posteriori error indicator, based on the stress recovery technique.

The layout of this paper is as follow: Section 2 provides a brief description of the constrained natural neighbor Galerkin method, applied to the elastostatic formulation. In Section 3, we will propose an a posteriori error indicator, as well as the associated refinement technique. Section 4 illustrates the capabilities of the proposed strategy through some numerical examples including different levels of stress concentration, all of them with known analytical solutions. Effectivity indexes will be computed to evaluate the accuracy of the error indicator.

## 2. The constrained natural element method

The natural element method is an appealing choice among the different meshless methods, because essential boundary conditions can be imposed directly, without determination to the other properties (linear consistency, smoothness, ...). It was originally proposed by Traversoni [39], Sambridge et al. [34] and was widely investigated in elastostatic problems by Sukumar [36]. Moreover, it was successfully applied in fluids dynamics simulations using updated Lagrangian formulations involving polymer injection and metal forming simulations [31,32,1], in problems involving cracks [42] and in phase change problems involving moving discontinuities [43].

### 2.1. Constrained natural element interpolation

We briefly touch upon the foundation of Sibson's natural neighbor coordinates (shape functions) that are used in the natural element method. For a more in-depth discussion on the Sibson interpolant and its application for solving second-order partial differential equations, the interested reader can refer to Braun and Sambridge [34] and Sukumar et al. [36]. The NEM interpolant is constructed on the underlying Voronoi diagram. The Delaunay tesselation is the topological dual of the Voronoi diagram.

For the sake of simplicity we only consider in the present work the 2D case, the 3D case being a direct extension.

Let $S = \{n_1, n_2, ..., n_N\}$ be a set of nodes in $\mathbb{R}^2$. The Voronoi diagram is the subdivision of $\mathbb{R}^2$ into regions $T_i$ (Voronoi cells) defined by

$$T_i = \{ \mathbf{x} \in \mathbb{R}^2 : d(\mathbf{x}, \mathbf{x}_i) < d(\mathbf{x}, \mathbf{x}_j), \forall j \neq i \}, \quad \forall i \quad (1)$$

The Sibson coordinates of $\mathbf{x}$ with respect to a natural neighbor $n_i$ (see Fig. 1) is defined as the ratio of the overlap area (volume in 3D) of their Voronoi cells to the total area

<figure><img src="image_1.png" /><figcaption>Fig. 1. Construction of the Sibson shape functions.</figcaption></figure>

(volume in 3D) of the Voronoi cell associated with the point $\mathbf{x}$:

$$\phi_i(\mathbf{x}) = \frac{\text{Area}(\text{afghe})}{\text{Area}(\text{abcde})} \quad (2)$$

If the point $\mathbf{x}$ coincides with the node $n_i$, i.e. $\mathbf{x} = \mathbf{x}_i$ then $\phi_i(\mathbf{x}_i) = 1$, and all the other shape functions vanish, i.e. $\phi_i(\mathbf{x}_j) = \delta_{ij}$ ($\delta_{ij}$ being the Kronecker's delta). The properties of positivity, interpolation and partition of unity are then verified [36]:

$$\left\{ \begin{array}{l} 0 \le \phi_i(\mathbf{x}) \le 1 \\ \phi_i(\mathbf{x}_j) = \delta_{ij} \\ \sum_{i=1}^{n} \phi_i(\mathbf{x}) = 1 \end{array} \right. \quad (3)$$

The natural neighbor interpolation satisfies the local coordinate property [35], namely:

$$\mathbf{x} = \sum_{i=1}^{n} \phi_i(\mathbf{x})\mathbf{x}_i \quad (4)$$

which combined with Eqs. (3), implies that the natural neighbor interpolant spans the space of linear polynomials (linear completeness).

Natural neighbor shape functions are $C^\infty$ at any point except at the nodes, where they are only $C^0$, and on the boundaries of the Delaunay circles (spheres in 3D) where they are only $C^1$, because of the discontinuity in the neighbor nodes across these boundaries. Hiyoshi and Sugihara [18] have shown that the Sibson interpolant belongs to a more general class of Voronoi-based interpolants, called $k$th order standard coordinates, proving that the interpolant generated by the $k$th order standard coordinates have $C^k$ continuity on the Delaunay circles (spheres) boundaries. Thus, the Sibsonian and non-Sibsonian (Laplace) coordinates [37] results to be the standard coordinates of order 1 and 0, respectively.

Another important property of this interpolation is its ability to reproduce linear functions on the boundary of convex domains. The proof can be found in Sukumar et al. [36], that we illustrate in Fig. 1 (b): due to the fact that the Voronoi cells are associated to points on the boundary become infinite, the contribution of internal points vanishes in the limit when the point approaches the convex boundary, and the shape functions associated with nodes $n_1$ and $n_2$ become linear on the segment $(n_1 - n_2)$. This is
not the case when non-convex boundaries are considered. This drawback will be considered later.

Consider an interpolation scheme for a vector-valued function $\mathbf{u}(\mathbf{x}) : \Omega \subset \mathbb{R}^2 \to \mathbb{R}^2$, in the form:

$$ \mathbf{u}^h(\mathbf{x}) = \sum_{i=1}^{n} \phi_i(\mathbf{x})\mathbf{u}_i \quad (5) $$

where $\mathbf{u}_i$ are the vectors defining the nodal degrees of freedom at the $n$ natural neighbors of point $\mathbf{x}$, and $\phi_i(\mathbf{x})$ are the Sibson coordinates defined in Eq. (2) associated with each node $n_i$. It can be noticed that Eq. (5) defines a local interpolation scheme which will be used to define both the trial and the test functions considered in the discretization of different variational formulations.

A recent development in the NEM, the constrained natural element method (C-NEM), was proposed in [42,43] in order to circumvent the problems induced by non-convex domains. In this approach, a visibility criterion is introduced to restrict the natural neighbors (influent nodes). For this purpose the constrained Voronoi diagram is constructed, from which the shape functions can be easily computed. In this manner, linear interpolation is recovered along the boundary of non-convex domains, making possible the introduction of essential boundary conditions as well as the treatment of fixed or moving discontinuities.

In the C-NEM framework, the interpolation can be expressed by

$$ \mathbf{u}^h(\mathbf{x}) = \sum_{i=1}^{V} \phi_i^C(\mathbf{x})\mathbf{u}_i \quad (6) $$

where $V$ is the number of natural neighbors visible from point $\mathbf{x}$ and $\phi_i^C$ is the constrained natural neighbor shape function, which is actually the Sibson interpolant computed using the constrained Voronoi diagram [42].

## 2.2. Elastostatic formulation

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

$$ \nabla \cdot \boldsymbol{\sigma} + \mathbf{b} = 0 \quad \text{in } \Omega \quad (7) $$

where $\Omega \in \mathbb{R}^2$ is the material domain, $\nabla$ is the divergence operator, $\boldsymbol{\sigma}$ is the Cauchy stress tensor and $\mathbf{b}$ is a body force term.

The constitutive relation is given by

$$ \boldsymbol{\sigma} = \mathbf{C}\boldsymbol{\epsilon} \quad (8) $$

where $\boldsymbol{\sigma}$ and $\boldsymbol{\epsilon}$ are the vector form of the stress and linearized strain tensor (symmetric part of the gradient of displacements tensor) respectively, and $\mathbf{C}$ is the elastic tensor matrix form.

The essential and natural boundary conditions are given by

$$ \begin{aligned} \mathbf{u} &= \bar{\mathbf{u}} & \text{on } \Gamma_u \\ \boldsymbol{\sigma}\mathbf{n} &= \bar{\mathbf{t}} & \text{on } \Gamma_t \end{aligned} \quad (9) $$

where $\Gamma = \Gamma_u \cup \Gamma_t$ is the boundary of $\Omega$, $\mathbf{n}$ is the unit outward normal vector defined on $\Gamma$, and $\bar{\mathbf{u}}$ and $\bar{\mathbf{t}}$ are the prescribed displacements and tractions, respectively.

The variational formulation associated with the elastostatic problem results:

Find $\mathbf{u} \in H^1(\Omega)$ kinematically admissible ($\mathbf{u} = \bar{\mathbf{u}}$ on $\Gamma_u$) such that

$$ \int_{\Omega} \boldsymbol{\sigma} \cdot \boldsymbol{\epsilon}^* \, d\Omega = \int_{\Omega} \mathbf{b} \cdot \mathbf{v}^* \, d\Omega + \int_{\Gamma_t} \bar{\mathbf{t}} \cdot \mathbf{v}^* \, d\Gamma, \quad \forall \mathbf{v}^* \in H_0^1(\Omega) \quad (10) $$

where $H^1(\Omega)$ and $H_0^1(\Omega)$ 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 $\mathbf{v}^*$, the following linear system of equations is obtained after numerical integration (whose details will be presented in the next section):

$$ \mathbf{K}\mathbf{d} = \mathbf{f}^{\text{ext}} \quad (11) $$

where $\mathbf{d}$ is the vector containing the nodal displacements, the matrix $\mathbf{K}$ and the vector $\mathbf{f}^{\text{ext}}$ being given by

$$ \mathbf{K} = \int_{\Omega} \mathbf{B}^{\text{T}} \mathbf{C} \mathbf{B} \, d\Omega \quad (12) $$

$$ \mathbf{f}^{\text{ext}} = \int_{\Gamma_t} \mathbf{N}^{\text{T}} \bar{\mathbf{t}} \, d\Gamma + \int_{\Omega} \mathbf{N}^{\text{T}} \mathbf{b} \, d\Omega \quad (13) $$

where $\mathbf{N}$ and $\mathbf{B}$ are the matrix containing the shape functions and their derivatives respectively.

## 2.3. Stabilized conforming numerical integration

Recently, new numerical integration procedures have been proposed in the framework of meshless simulations. Domain integration using Gauss quadrature introduces significant numerical errors due to the following statements: (i) the integration cells do not match shape function supports; (ii) the non-polynomial character of the meshless shape functions [13]. On the other hand, direct nodal integration, using the nodes as integration points, leads to numerical instabilities [9].

In this work, we use the integration procedure proposed by Chen et al. [8]: the stabilized conforming nodal integration (SCNI), which was applied to the NEM by Gonzalez et al. [15] and Yoo et al. [40]. In this approach, a strain smoothing stabilization is performed to stabilize the nodal integration. The SCNI is based on the assumed strain method, in which a modified gradient is introduced at the integration point (node) [8]:

$$ \tilde{\mathbf{V}}\mathbf{u}^h(\mathbf{x}_i) = \frac{1}{A_i} \int_{\Omega_i} \mathbf{V}\mathbf{u}^h(\mathbf{x}) \, d\Omega \quad (14) $$

where $\mathbf{x}_i$ are the coordinates of node $n_i$.
Thus the modified strain vector is given by

$$ \tilde{\epsilon}^h(\mathbf{x}_i) = \frac{1}{A_i} \int_{\Omega_i} \epsilon^h(\mathbf{x}) \, d\Omega = \frac{1}{A_i} \int_{\Omega_i} \left\{ \begin{matrix} \frac{\partial u_1^h(\mathbf{x})}{\partial x_1} \\ \frac{\partial u_2^h(\mathbf{x})}{\partial x_2} \\ \frac{\partial u_1^h(\mathbf{x})}{\partial x_2} + \frac{\partial u_2^h(\mathbf{x})}{\partial x_1} \end{matrix} \right\} d\Omega \quad (15) $$

Now, by applying the divergence theorem, it results:

$$ \tilde{\epsilon}^h(\mathbf{x}_i) = \frac{1}{A_i} \int_{\Gamma_i} \left\{ \begin{matrix} u_1^h(\mathbf{x})n_1 \\ u_2^h(\mathbf{x})n_2 \\ u_1^h(\mathbf{x})n_2 + u_2^h(\mathbf{x})n_1 \end{matrix} \right\} d\Gamma \quad (16) $$

$\Omega_i$ is the Voronoi cell related to the node $n_i$ whose boundary will be noted by $\Gamma_i$, being $A_i$ its area, as depicted in Fig. 2.

Introducing now the natural element approximation of the displacement field according to Eq. (6) we obtain:

$$ \tilde{\epsilon}^h(\mathbf{x}_i) = \tilde{\mathbf{B}}_i d \quad (17) $$

which can be expressed explicitly by

$$ \tilde{\epsilon}^h(\mathbf{x}_i) = \begin{pmatrix} \frac{\partial \tilde{\phi}_1(\mathbf{x}_i)}{\partial x_1} & 0 & \frac{\partial \tilde{\phi}_2(\mathbf{x}_i)}{\partial x_1} & 0 & \dots & \frac{\partial \tilde{\phi}_N(\mathbf{x}_i)}{\partial x_1} & 0 \\ 0 & \frac{\partial \tilde{\phi}_1(\mathbf{x}_i)}{\partial x_2} & 0 & \frac{\partial \tilde{\phi}_2(\mathbf{x}_i)}{\partial x_2} & \dots & 0 & \frac{\partial \tilde{\phi}_N(\mathbf{x}_i)}{\partial x_2} \\ \frac{\partial \tilde{\phi}_1(\mathbf{x}_i)}{\partial x_2} & \frac{\partial \tilde{\phi}_1(\mathbf{x}_i)}{\partial x_1} & \frac{\partial \tilde{\phi}_2(\mathbf{x}_i)}{\partial x_2} & \frac{\partial \tilde{\phi}_2(\mathbf{x}_i)}{\partial x_1} & \dots & \frac{\partial \tilde{\phi}_N(\mathbf{x}_i)}{\partial x_1} & \frac{\partial \tilde{\phi}_N(\mathbf{x}_i)}{\partial x_2} \end{pmatrix} = \tilde{\mathbf{B}}_i d \quad (18) $$

where the components of $\tilde{\mathbf{B}}_i$ are defined by

$$ \frac{\partial \tilde{\phi}_j(\mathbf{x}_i)}{\partial x_1} = \frac{1}{A_i} \int_{\Gamma_i} \phi_j(\mathbf{x}) n_1(\mathbf{x}) d\Gamma \quad (19) $$

$$ \frac{\partial \tilde{\phi}_j(\mathbf{x}_i)}{\partial x_2} = \frac{1}{A_i} \int_{\Gamma_i} \phi_j(\mathbf{x}) n_2(\mathbf{x}) d\Gamma \quad (20) $$

Obviously, most of the components in the matrix $\tilde{\mathbf{B}}_i$ are zero due to the compact support of the shape functions. Chen et al. [8] proposed to use the intersection between the Voronoi diagram and the domain closure to define the integration cells $\Omega_i$. 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_i$:

$$ \mathbf{K} = \sum_{i} \tilde{\mathbf{K}}_{i} = \sum_{i} A_{i} \tilde{\mathbf{B}}_{i}^{\mathrm{T}} \tilde{\mathbf{C}} \tilde{\mathbf{B}}_{i} \quad (21) $$

and then, the stress in each cell can be deduced from:

$$ \tilde{\sigma}_{i}^{h} = \tilde{\mathbf{C}} \tilde{\mathbf{B}}_{i} d \quad (22) $$

which results also constant inside each Voronoi cell.

The use of the divergence theorem avoids the computation of the shape functions derivatives, requiring only the evaluation of the shape functions on the boundary of the Voronoi cells.

<figure><img src="image_2.png" /><figcaption>Fig. 2. Integration domains for the SCNI procedure: (a) intersection between the Voronoi diagram and the domain closure and (b) constrained and bounded Voronoi diagram.</figcaption></figure>

This integration technique can be naturally applied in the C-NEM context, where the nodal integration domains coincide with the intersection between the constrained Voronoi cells and the domain. Moreover, it must be noticed, as illustrated in Fig. 2, that the intersection between the standard Voronoi diagram and the domain does not define accurately the integration domain.

### 3. Adaptation procedure

Error indicators are used in numerical simulations to provide some information about the refinement or nodal adaptation to carry out, for improving the solution accuracy. In this section, an error indicator is proposed and tested in the C-NEM context, that will be used to define an adaptation strategy.

#### 3.1. Error indicator based on natural neighbor recovery stress

The Zienkiewicz–Zhu estimator was proposed in [47]. Its aim is to estimate $||\mathbf{V}\mathbf{u} - \mathbf{V}\mathbf{u}^h||_{L^2(\Omega)}$. For this purpose, a higher-order recovery $\mathbf{G}\mathbf{u}$ of $\mathbf{V}\mathbf{u}$ must be constructed using only $\mathbf{u}^h$. We first consider an approximation $\mathbf{G}\mathbf{u}$ more accurate than $\mathbf{V}\mathbf{u}^h$, i.e.

$$ \|\mathbf{V}\mathbf{u} - \mathbf{G}\mathbf{u}\|_{L^2(\Omega)} \le C_1 h^\alpha \|\mathbf{V}\mathbf{u} - \mathbf{V}\mathbf{u}^h\|_{L^2(\Omega)}, \quad \alpha > 0 \quad (23) $$

The effectivity index is a measure of the error estimate compared to the exact error. We have

$$ \|\mathbf{V}\mathbf{u}^h - \mathbf{G}\mathbf{u}\|_{L^2(\Omega)} = \|(\mathbf{V}\mathbf{u}^h - \mathbf{V}\mathbf{u}) - (\mathbf{G}\mathbf{u} - \mathbf{V}\mathbf{u})\|_{L^2(\Omega)} \quad (24) $$

Using triangle inequality, we have

$$ \begin{aligned} & \|\mathbf{V}\mathbf{u}^h - \mathbf{V}\mathbf{u}\|_{L^2(\Omega)} - \|\mathbf{G}\mathbf{u} - \mathbf{V}\mathbf{u}\|_{L^2(\Omega)} \le \|\mathbf{V}\mathbf{u}^h - \mathbf{G}\mathbf{u}\|_{L^2(\Omega)} \\ & \le \|\mathbf{V}\mathbf{u}^h - \mathbf{V}\mathbf{u}\|_{L^2(\Omega)} + \|\mathbf{G}\mathbf{u} - \mathbf{V}\mathbf{u}\|_{L^2(\Omega)} \end{aligned} \quad (25) $$
Then dividing each term by $||\mathbf{V}\mathbf{u}^h - \mathbf{V}\mathbf{u}||_{L^2(\Omega)}$ it leads

$$1 - \frac{||\mathbf{G}\mathbf{u} - \mathbf{V}\mathbf{u}||_{L^2(\Omega)}}{||\mathbf{V}\mathbf{u}^h - \mathbf{V}\mathbf{u}||_{L^2(\Omega)}} \le \theta \le 1 + \frac{||\mathbf{G}\mathbf{u} - \mathbf{V}\mathbf{u}||_{L^2(\Omega)}}{||\mathbf{V}\mathbf{u}^h - \mathbf{V}\mathbf{u}||_{L^2(\Omega)}} \quad (26)$$

If Eq. (23) is verified, then the effectivity index is supposed approaching 1 as $h$ approaches zero, and thus leads to asymptotically exact error indication. Let $\psi_i(\mathbf{x})$ and $\phi_i(\mathbf{x})$ two sets of shape functions with zero and first order of consistency, respectively. We have for $\alpha > 0$:

$$\left\| \mathbf{V}\mathbf{u} - \sum_{i=1}^{n} \phi_i(\mathbf{x}) \mathbf{V}\mathbf{u}^h(\mathbf{x}_i) \right\|_{L^2(\Omega)} \le C_1 h^\alpha \left\| \mathbf{V}\mathbf{u} - \sum_{i=1}^{n} \psi_i(\mathbf{x}) \mathbf{V}\mathbf{u}^h(\mathbf{x}_i) \right\|_{L^2(\Omega)}$$

By choosing $\psi_i(\mathbf{x})$ as the Thiessen natural neighbor shape function

$$\psi_i(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \in \Omega_i \\ 0 & \text{otherwise} \end{cases} \quad (27)$$

where $\Omega_i$ is a Voronoi cell. By choosing $\phi_i(\mathbf{x})$ as the Sibson natural neighbor shape functions and by using the stabilized conforming nodal integration scheme previously described, we define

$$\tilde{\mathbf{V}}\mathbf{u}^h(\mathbf{x}_i) \approx \frac{1}{\Omega_i} \int_{\Omega_i} \mathbf{V}\mathbf{u}^h(\mathbf{x}) \, d\Omega \quad (28)$$

$$\tilde{\mathbf{V}}\mathbf{u}^h(\mathbf{x}) = \sum_{i=1}^{n} \psi_i(\mathbf{x}) \tilde{\mathbf{V}}\mathbf{u}^h(\mathbf{x}_i) \quad (29)$$

We then construct **Gu** using:

$$\mathbf{G}\mathbf{u}(\mathbf{x}) = \sum_{i=1}^{n} \phi_i(\mathbf{x}) \tilde{\mathbf{V}}\mathbf{u}^h(\mathbf{x}_i) \quad (30)$$

Using Eqs. (30), (29) and (23) we have

$$||\mathbf{V}\mathbf{u} - \mathbf{G}\mathbf{u}||_{L^2(\Omega)} \le C_1 h^\alpha ||\mathbf{V}\mathbf{u} - \tilde{\mathbf{V}}\mathbf{u}^h||_{L^2(\Omega)}, \quad \alpha > 0 \quad (31)$$

Now we define

$$\tilde{\sigma}^h(\mathbf{x}) = C\tilde{\mathbf{B}}d \quad (32)$$

with $\tilde{\mathbf{B}}$ defined in Eq. (18), as well as

$$\hat{\sigma}(\mathbf{x}) = \sum_{i=1}^{n} \phi_i(\mathbf{x}) \tilde{\sigma}^h(\mathbf{x}_i) \quad (33)$$

In the case of linear problems, $\tilde{\mathbf{B}}$ is linear with respect to $\tilde{\sigma}^h(\mathbf{x})$. We thus have

$$||\mathbf{u}^{\tilde{\sigma}}(\mathbf{x}) - \mathbf{u}^{\sigma}(\mathbf{x})||_{E(\Omega)} \le C_1 h^{\alpha} ||\mathbf{u}^{\tilde{\sigma}}(\mathbf{x}) - \mathbf{u}^{\sigma}(\mathbf{x})||_{E(\Omega)}, \quad \alpha > 0 \quad (34)$$

where $||\mathbf{u}^{\sigma}(\mathbf{x})||_{E(\Omega)}$ denotes the energy norm with respect to $\mathbf{u}(\mathbf{x})$, defined as

$$||\mathbf{u}^{\sigma}(\mathbf{x})||_{E(\Omega)} = \left\{ \int_{\Omega} \frac{1}{2} \boldsymbol{\sigma}^{\mathrm{T}} \mathbf{C}^{-1} \boldsymbol{\sigma} \right\}^{1/2} \quad (35)$$

<figure><img src="image_3.png" /><figcaption>Fig. 3. Computing errors in each Voronoi cell.</figcaption></figure>

Finally, using (26) and (34), we have

$$1 - \frac{||\mathbf{u}^{\tilde{\sigma}} - \mathbf{u}^{\sigma}||_{E(\Omega)}}{||\mathbf{u}^{\tilde{\sigma}} - \mathbf{u}^{\sigma}||_{E(\Omega)}} \le \theta \le 1 + \frac{||\mathbf{u}^{\tilde{\sigma}} - \mathbf{u}^{\sigma}||_{E(\Omega)}}{||\mathbf{u}^{\tilde{\sigma}} - \mathbf{u}^{\sigma}||_{E(\Omega)}} \quad (36)$$

Then $\theta$ is an effectivity index associated with the proposed error indicator, assumed to be asymptotically exact.

**Remark.** The choice of $\phi_i(\mathbf{x})$ is not unique, one could consider higher-order reproducing shape functions. For example, we have proposed in [45] a quadratic NEM interpolant, that may lead to more accurate error indicator. Nevertheless, the higher cost of that interpolant may not be worth as compared to the simplicity of the proposed indicator. A study of such indicator is outside the scope of this paper.

In order to evaluate the integral in Eq. (35), the Voronoi cells are triangulated, as depicted in Fig. 3, and Gauss integration is performed in each triangle. As the constrained Voronoi cells define a partition of the whole domain, the global error can be expressed as a function of local contributions by

$$||\mathbf{u}^{\tilde{\sigma}} - \mathbf{u}^{\sigma}||_{E(\Omega)}^2 = \sum_{i=1}^{N} ||\mathbf{u}^{\tilde{\sigma}} - \mathbf{u}^{\sigma}||_{E(\Omega_i)}^2 \quad (37)$$

where $N$ is the total number of nodes in the domain. We can also define a local relative contribution to the global error from:

$$\eta_i = \frac{||\mathbf{u}^{\tilde{\sigma}} - \mathbf{u}^{\sigma}||_{E(\Omega_i)}}{||\mathbf{u}^{\tilde{\sigma}}||_{E(\Omega)}} \quad (38)$$

### 3.2. Refinement strategy

To improve the solution approximation without a detriment in the computational efficiency, meanwhile to adapt the nodal density and distribution according to the solution behavior, the knowledge of local error contributions is crucial.

For this purpose, we define a permissible error indication $\tilde{\eta}$ and the permissible error indication in each cell $\tilde{\eta}_i$ satisfying:
$$ \bar{\eta}^2 = \sum_{i=1}^{N} (\bar{\eta}_i)^2 \qquad (39) $$

If we are looking for a similar amount of error in each cell, then $\bar{\eta}_i = \bar{\eta}_*$, $\forall i$:

$$ \bar{\eta}_* = \frac{\bar{\eta}}{\sqrt{N}} \qquad (40) $$

Following Zienckiewicz et al. [48], we can compute the characteristic length of the new cell with respect to the old one by assuming a p-order convergence ($\Theta(h^p)$). Thus, we can write

$$ \frac{h_i^{\text{new}}}{h_i^{\text{old}}} = \left[ \frac{\bar{\eta}_*}{\bar{\eta}_i} \right]^{\frac{1}{p}} = \left[ \frac{\bar{\eta} \|\mathbf{u}^\sigma\|_{E(\Omega)}}{\sqrt{N} \|\mathbf{u}^\sigma - \mathbf{u}^\sigma\|_{E(\Omega_i)}} \right]^{\frac{1}{p}} \qquad (41) $$

From our experience, the C-NEM rate of convergence is very close to that obtained using triangular finite elements with linear interpolation (due to the linear consistency of Sibson's interpolation). Thus, p is expected to be the rate of convergence using an energy norm.

To decide whether a Voronoi cell should be subdivided, we define the characteristic length h of a Voronoi cell as the mean distance between the central node and its visible neighbors, which using any of the refinement strategies shown in Fig. 4, leads to

$$ \frac{h_i^{\text{new}}}{h_i^{\text{old}}} \le \frac{1}{2} \qquad (42) $$

Thus, the refinement must operate if

$$ \left[ \frac{\bar{\eta} \|\mathbf{u}^\sigma\|_{E(\Omega)}}{\sqrt{N} \|\mathbf{u}^\sigma - \mathbf{u}^\sigma\|_{E(\Omega_i)}} \right]^{\frac{1}{p}} \le \frac{1}{2} \qquad (43) $$

The strategy depicted in Fig. 4(a) has been considered by several authors [41,29] and is suitable to refine irregular clouds of nodes. Starting from an initial cloud of nodes, we first construct the associated Voronoi diagram. If the condition (43) is valid in a given Voronoi cell, then new nodes are added on the vertices of the cell. This operation is repeated for all cells. Then We reconstruct the Voronoi diagram of the new cloud of nodes, and so on until the global error is below the given tolerance.

<figure><img src="image_4.png" /><figcaption>Fig. 4. Strategies of nodal refinement based on the Voronoi cells. (o): additional nodes.</figcaption></figure>

A different strategy is described in Fig. 4(b), where the new nodes are inserted between the central node and the natural (visible) neighbors. In the following, this strategy (b) has been applied to treat the numerical examples.

## 4. Numerical examples

### 4.1. Plane stress plate problem

A square plate of length L is subjected to the following system of loads (as depicted in Fig. 5):

$$ F_d = \left( \frac{2x_1}{L} - 1 \right) e_1 - \frac{x_2}{L} \left( 1 - \frac{x_2}{L} \right) e_2 \quad \text{on } x_1 = 0 \qquad (44) $$

and

<figure><img src="image_5.png" /><figcaption>Fig. 5. (a) Square plate with plane stresses. (b) Constrained Voronoi cells.</figcaption></figure>
$$F_d = \frac{x_2}{L} \left(1 - \frac{x_2}{L}\right) e_2 \quad \text{on } x_1 = L \qquad (45)$$

where $e_1$ and $e_2$ are the unit vectors defining the $x_1$ and the $x_2$ axes directions.

This benchmark was investigated in the finite element context by Ladevèze and Pelle [23]. The exact stress solution is given by

$$\sigma_{11} = \left(\frac{2x_2}{L} - 1\right) \left(\frac{x_1}{L} - 1\right) \qquad (46)$$

$$\sigma_{22} = 0 \qquad (47)$$

$$\sigma_{12} = \frac{x_2}{L} \left(1 - \frac{x_2}{L}\right) \qquad (48)$$

In order to evaluate the error indicators, the energy norms related to the exact stress and the different computed stresses are plotted in Fig. 6 for different nodal densities. We can notice that the reference solution $\hat{\sigma}$ converges faster than $\tilde{\sigma}$. The effectivity index is then expected to approach asymptotically to one in this particular problem. The energy norm of the different errors indications are depicted in Fig. 7, from which we can notice the higher rate of convergence of $||\mathbf{u}^{\hat{\sigma}} - \mathbf{u}^{\sigma^{\text{ex}}}||_{E(\Omega)}$.

Fig. 8 and Table 1 present the effectivity index associated with the proposed error indicator. The index is close to one, and the accuracy increases asymptotically, as expected. Nevertheless, the stress field is smooth in this problem, and in order to evaluate the error indicator for nodal refinement purpose, we will study in the next section problems involving stress concentration.

## 4.2. Traction of a plate with a central hole

A traction load is applied on a square plate as depicted in Fig. 9 under plane stress conditions. The exact stress solution for an infinite plate subjected to a unit traction on the boundaries whose unit normal vector coincides with the $x_1$ direction, is given by [38]

<figure><img src="image_6.png" /><figcaption>Fig. 6. Convergence of the exact and computed stresses for the plane stress plate problem.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7. Convergence of the different errors in a plane stress plate problem.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 8. Effectivity index associated with the NN1 and the NN2 error indicators for the plane stress plate problem.</figcaption></figure>

<table><caption>Table 1
Effectivity indexes for the different problems</caption><thead><tr><th rowspan="2">Nb. nodes</th><th>Example 1</th><th>Example 2</th><th>Example 3</th></tr><tr><th>Smooth stress field</th><th>Stress concentration</th><th>Stress singularity</th></tr></thead><tbody><tr><td>5 × 5</td><td>0.925</td><td>0.746</td><td>0.721</td></tr><tr><td>10 × 10</td><td>0.967</td><td>0.891</td><td>0.748</td></tr><tr><td>20 × 20</td><td>0.984</td><td>0.969</td><td>0.722</td></tr><tr><td>30 × 30</td><td>0.989</td><td>0.969</td><td>0.720</td></tr></tbody></table>

$$\sigma_{ex11} = 1 - \frac{a^2}{r^2} \left\{ \frac{3}{2} \cos(2\theta) + \cos(4\theta) \right\} + \frac{3}{2} \frac{a^4}{r^4} \cos(4\theta) \quad (49)$$

$$\sigma_{ex22} = -\frac{a^2}{r^2} \left\{ \frac{1}{2} \cos(2\theta) - \cos(4\theta) \right\} - \frac{3}{2} \frac{a^4}{r^4} \cos(4\theta) \quad (50)$$

$$\sigma_{ex12} = -\frac{a^2}{r^2} \left\{ \frac{1}{2} \sin(2\theta) + \sin(4\theta) \right\} + \frac{3}{2} \frac{a^4}{r^4} \sin(4\theta) \quad (51)$$
<figure><img src="image_9.png" /><figcaption>Fig. 9. (a) Square plate with a central hole. (b) Computational domain and boundary conditions.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 10. Energy norm of the errors in the plate with a hole problem.</figcaption></figure>

To avoid the error introduced by the finite size of the plate, exact traction $F_d = \sigma_{ex}n$ are imposed on the square plate boundaries. For reasons of symmetry, only a quarter of the plate is analyzed. From Fig. 10 we can notice that the reference solution $\hat{\sigma}$ still converges faster than other ones. The effectivity index is presented in Fig. 11 as well as in Table 1, exhibiting the asymptotically increase of accuracy, despite of the stress concentration. A refinement sequence using the criterion (43) is depicted in Fig. 12. As we can notice the refinement is activated in the regions with a significant stress concentration.

### 4.3. Mode-I crack problem

In this last example, a mode-I crack problem is considered, whose geometry is depicted in Fig. 13(a).

The exact stress solution in the crack tip neighborhood is given by [26]

<figure><img src="image_11.png" /><figcaption>Fig. 11. Effectivity index in the plate with a hole problem.</figcaption></figure>

$$ \sigma_{11} = \frac{K_I}{\sqrt{2\pi r}} \cos\frac{\theta}{2} \left( 1 - \sin\frac{\theta}{2} \sin\frac{3\theta}{2} \right) \quad (52) $$

$$ \sigma_{22} = \frac{K_I}{\sqrt{2\pi r}} \cos\frac{\theta}{2} \left( 1 + \sin\frac{\theta}{2} \sin\frac{3\theta}{2} \right) \quad (53) $$

$$ \sigma_{12} = \frac{K_I}{\sqrt{2\pi r}} \sin\frac{\theta}{2} \cos\frac{\theta}{2} \cos\frac{3\theta}{2} \quad (54) $$

To extend this solution to the whole domain, we impose on the square plate boundary the exact traction: $F_d = \sigma_{ex}n$, where we assume $L = 1$ mm, $a = L/2$ mm and $K_I = 1$. The constrained Voronoi diagram allow us to introduce easily and accurately the presence of the crack as well as continuity around the crack tip without any special treatment due to the salient features of the C-NEM shape functions support [42].

The effectivity index is plotted in Fig. 14. We can notice that the error indicator is reasonable and remain constant for different nodal densities, perhaps due to the fact that the stabilized conforming nodal integration regularizes
<figure><img src="image_12.png" /><figcaption>Fig. 12. Adaptive refinement in the plate with a hole problem.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 13. (a) Mode-I crack problem geometry. (b) Exact boundary conditions.</figcaption></figure>

the singularity, making it a useful tool for nodal refinement.

A sequence of the adaptive refinement using the proposed error indicator is depicted in Fig. 15. We can observe
<figure><img src="image_14.png" /><figcaption>Fig. 14. Effectivity index in the mode-I crack problem.</figcaption></figure>

that the stress concentration in the singularity neighborhood is accurately accounted.

## 5. Conclusions

In this paper a refinement strategy able to operate with accuracy in the C-NEM framework was presented and has been tested in some linear elastostatic problems. A simple a posteriori error indicator based on recovery stress in tandem with stabilized nodal integration has been proposed and tested in some benchmarks problems involving different levels of difficulty (stress concentration, solution singularity, ...) all of them with a known exact solution. The indicator exhibits asymptotically exact error estimation in problems with different level of stress concentrations. In problems involving singularities, the effectivity is reasonable and constant, which makes it a useful tool for nodal refinement in the context of the natural element method.

<figure><img src="image_15.png" /><figcaption>Fig. 15. Refinement procedure in the mode-I crack problem.</figcaption></figure>
Different adaptation procedures have been performed successfully, proving the utility of the procedure.

## References

[1] Alfaro I, Yvonnet J, Cueto E, Doblaré M, Chinesta F. Meshless methods with application to metal forming. Comput Methods Appl Mech Eng, in press.

[2] Babuška I, Aziz AK. On the angle condition in the finite element method. SIAM J Numer Anal 1976;13(2):214-26.

[3] Babuška I, Rheinboldt WC. Adaptive approaches and reliability estimations in finite element analysis. Comput Methods Appl Mech Eng 1979;17/18:519-40.

[4] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37:229-56.

[5] Belytschko T, Liu WK, Singer M. On adaptivity and error criteria for meshfree methods. In: Ladevèze P, Oden JT, editors. Advances in adaptive computational methods in mechanics, 1998. p. 217-28.

[6] Chen JS, Wang HP. New boundary conditions treatments in meshfree computation of contact problems. Comput Methods Appl Mech Eng 2000;187:441-68.

[7] Chen JS, Wu CT, Yoon Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. Int J Numer Methods Eng 2001;50: 435-66.

[8] Chen JS, Yoon S, Wu CT. Non-linear version of stabilized conforming nodal integration for Galerkin mesh-free methods. Int J Numer Methods Eng 2002;53:2587-615.

[9] Chung H-J, Belytschko T. An error estimate in the EFG method. Comput Mech 1998;21:91-100.

[10] Cueto E, Cegoñino J, Calvo B, Doblaré M. On the imposition of essential boundary conditions in natural neighbor Galerkin methods. Commun Numer Methods Eng 2003;19:361-76.

[11] De S, Bathe KJ. The method of finite spheres with improved numerical integration. Comput Struct 2005;79:2183-96.

[12] Dolbow J, Belytschko T. Numerical integration of the Galerkin weak form in meshfree methods. Comput Mech 1999;23:219-30.

[13] Duarte CA, Oden JT. An H-p adaptative method using clouds. Comput Methods Appl Mech Eng 1996;139:237-62.

[14] Gonzalez D, Cueto E, Martinez MA, Doblaré M. Numerical integration in natural neighbor Galerkin methods. Int J Numer Methods Eng 2002;60(12):2077-104.

[15] Grätsch T, Bathe KJ. A posteriori error estimation techniques in practical finite element analysis. Comput Struct 2005;83:235-65.

[16] Hiyoshi H, Sugura K. Improving continuity of Voronoi-based interpolation over Delaunay spheres. Comput Geometry 2002;22: 167-83.

[17] Hong JW, Bathe KJ. Coupling and enrichment schemes for finite element and finite sphere discretizations. Comput Struct 2005;83: 1386-95.

[18] Krok J, Orkisz J. Application of the generalized FDM to calculation of arbitrary loaded axisymmetrical massive structures. In: Proc of 28th conf. KILiW PAN and KN PZITB, Krynica, Poland, 1982. p. 81-90 [1992, vol. 10, p. 307-18].

[19] Ladevèze P, Pelle JP, Rougeot P. Error estimation and mesh optimization for classical finite elements. Eng Comput 1991;8:69-80.

[20] Ladevèze P, Pelle J-P. La Maitrise du calcul en mécanique linéaire et non-linéaire. Hermes/Lavoisier, 2001.

[21] Lee CK, Zhou CE. On error estimation and adaptive refinement for element free Galerkin method: Part I: stress recovery and a posteriori error estimation. Comput Struct 2003;82(4-5):413-28.

[22] Lee CK, Zhou CE. On error estimation and adaptive refinement for element free Galerkin method: Part II: adaptive refinement. Comput Struct 2003;82(4-5):429-43.

[23] Perrone N, Liebowitz H, Mulville D, Plikey W. Fracture mechanics. Charlottesville: University Press of Virginia; 1978.

[24] Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. Int J Numer Methods Fluids 1995;21:1081-106.

[25] Liu WK, Uras RA, Chen Y. Enrichment of the finite element method with the reproducing kernel particle method. J Appl Mech ASME 1997;64(4):861-70.

[26] Lu H, Chen JS. Adaptive meshfree particle method. Lect Notes Comput Sci Eng 2002;26:251-67.

[27] Lucy LB. A numerical approach to the testing of fusion process. Astronom J 1977;88:1013-24.

[28] Martinez MA, Cueto E, Doblaré M, Chinesta F. A meshless simulation of injection processes involving short fibers composites. Int J Forming Process 2001;4/3:217-36.

[29] Martinez MA, Cueto E, Doblaré M, Chinesta F. Updated lagrangian free surface flow simulations with natural neighbour Galerkin methods. Int J Numer Methods Eng 2004;60(13):2105-29.

[30] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992;10:307-18.

[31] Sambridge M, Braun J, McQueen M. Geophysical parameterization and interpolation of irregular data using natural neighbors. Geophys J Int 1995;122:837-57.

[32] Sibson R. A vector identity for the Dirichlet tesselations. Math Proc Camb Phil Soc 1980;87:151-5.

[33] Sukumar N, Moran B, Belytschko T. The natural elements method in solid mechanics. Int J Numer Methods Eng 1998;43: 839-87.

[34] Sukumar N, Moran B, Semenov Y, Belikov VV. Natural neighbor Galerkin methods. Int J Numer Methods Eng 2001;50:207-19.

[35] Timoshenko SP, Goodier JN. Theory of elasticity. 3rd ed. New York: McGraw Hill; 1970.

[36] Traversoni L. Natural neighbor finite elements. In: International conference on hydraulic engineering software. Hydrosoft Proceedings, vol. 2. Computational Mechanics Publications, p. 291-7.

[37] Yoon J, Moran B, Chen J-S. Stabilized conforming nodal integration in the natural-element method. Int J Numer Methods Eng 1998;60:861-90.

[38] You Y, Chen JS, Lu H. Filters, reproducing kernel, and adaptive meshfree method. Comput Mech 2003;31(3-4):316-26.

[39] Yvonnet J, Ryckelynck D, Lorong P, Chinesta F. A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM). Int J Numer Methods Eng 2004;60:1451-74.

[40] Yvonnet J, Chinesta F, Lorong P, Rynckelynck D. The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces. Int J Thermal Sci 2005;44: 559-69.

[41] Yvonnet J, Villon P, Chinesta F. Natural element approximations involving bubbles for treating mechanical models in incompressible media. Int J Numer Methods Eng 2006;66(7):1125-52.

[42] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. Int J Numer Methods Eng 1987;142:337-57.

[43] Zienkiewicz OC, Boroomand B, Zhu JZ. Recovery procedures in error estimation and adaptivity Part I: Adaptivity in linear problems. Comput Methods Appl Mech Eng 1999;176:111-25.