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

Julien Yvonnet, David Ryckelynck, Philippe Lorong, Francisco Chinesta

## ► To cite this version:

Julien Yvonnet, David Ryckelynck, Philippe Lorong, Francisco Chinesta. A new extension of the natural element method for non-convex and discontinuous problems: the constrained natural element method (C-NEM). International Journal for Numerical Methods in Engineering, 2004, 60 (8), pp.1451-1474. 10.1002/nme.1016. hal-01508695

HAL Id: hal-01508695

https://hal.science/hal-01508695v1

Submitted on 15 Apr 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

<figure><img src="image_1.png" /></figure>
# A new extension of the natural element method for non-convex and discontinuous problems: the constrained natural element method (C-NEM)

J. Yvonnet, D. Ryckelynck, P. Lorong and F. Chinesta

Laboratoire de Mécanique des Systèmes et des Procédés, UMR 8106 CNRS ENSAM-ESEM,
151 boulevard de l'Hôpital, Paris F-75013, France

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

KEY WORDS: C-NEM; non-convex domains; visibility criterion; constrained Voronoi diagram

## 1. INTRODUCTION

In recent years, a new family of computational methods have emerged. The so-called mesh-
less or mesh-free methods have been investigated and used by many researchers for treating
a large variety of engineering problems, involving usually large displacements as encountered
for example in forming process simulations (free surface or moving boundary problems, mov-
ing interfaces, cracks propagation, etc.). In these problems, accurate finite element solutions
require significant computational efforts in remeshing steps. In contrast, meshless methods re-
quire only nodal data without explicit connectivity between nodes. The first generation of
meshless methods, the smooth particle hydrodynamics (SPH) method, originally proposed by
Lucy [1], introduces kernel approximations and imposes the conservation laws using a collo-
cation technique. Recently, other meshfree methods have been proposed, such as the diffuse
element method (DEM) [2], the element free Galerkin (EFG) [3], the partition of unity method
(PUM) [4], the reproducing kernel particle method (RKPM) [5] or the HP-clouds method [6]
and they have been successfully used to discretize the Galerkin weak form associated with
several problems. Galerkin meshfree methods differ from SPH in two ways: (i) consistency is
introduced from kernel corrections, and (ii) discretization is applied to the weak formulation.

Nevertheless, a troublesome task in meshless methods is the imposition of essential boundary
conditions because the shape functions related to these methods do not satisfy the Kronecker
delta property and they do not vanish on the boundary where the essential conditions are
imposed. Special strategies are hence required to enforce the essential boundary conditions,
involving additional and costly computational efforts. Among them are the Lagrange multiplier
method [3], the transformation method [7], the approach based on the D'Alembert's principle
[8], the introduction of a singular weight function [9], the penalty method [10], the mixed
transformation method [11] and the boundary singular kernel method [11]. In order to avoid
this additional treatment, a new meshless method has recently been investigated: the natural
element method (NEM) [12–14]. In the NEM, the test and trial functions are constructed using
natural neighbour (Sibson) co-ordinates [15]. Thus, shape functions satisfy the Kronecker delta
property, allowing the imposition of the essential boundary conditions like in the finite element
method (FEM). Another important property of the natural neighbour interpolation is that it
satisfies linear consistency [15]. The natural neighbour co-ordinates are based on well-known
geometric concepts: the Voronoi diagram and the Delaunay tesselation. As the Voronoi diagram
and its dual Delaunay tesselation are only defined for convex domains, issues occur when non-
convex domains are considered. Thus, mutual influences between nodes on the boundary of
non-convex domains appear, which introduces difficulties in the imposition of essential boundary
conditions.

An extension of the NEM to facilitate its application for any domain (convex or not) and
retaining its attractive properties, the constrained natural element method (C-NEM) was proposed
in our former work [16]. The main originality of this approach is the introduction of a visibility
criterion [17] in the NEM to select natural neighbours for constructing the interpolation. A
modified, so-called constrained Voronoi diagram is introduced for the computation of the shape
functions. In addition, the stabilized conforming nodal integration (SCNI) scheme [18, 19] is
used instead of Gauss integration to improve both the accuracy and the efficiency of the method.

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

2. THE NATURAL ELEMENT METHOD

2.1. Voronoi diagram and natural neighbours

The notion of natural neighbours and natural neighbours interpolation were introduced by Sibson
[15] for data fitting and smoothing of highly irregular set of nodes. This interpolant is based
<figure><img src="image_2.png" /><figcaption>Figure 1. Voronoi diagram of a set of nodes in the plane.</figcaption></figure>

on geometric constructs known as Voronoi diagram and Delaunay tesselations. The concept of Voronoi diagram was originally introduced by mathematicians (Dirichlet, 1850; Voronoi, 1908) [20] and later applied in many fields of engineering sciences. In 2D, the Voronoi diagram of an irregular set of nodes divides the plane into a set of regions $T_i$, one for each node $n_i$, such that all points in a particular region are closer to its defining node than to any other node (see Figure 1). The Voronoi diagram is unique, it fills the convex hull of the set of nodes, and can be extended to any dimension. In mathematical terms, the Voronoi diagram is defined as

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

The other closely related geometric construct which is of interest here is the Delaunay tesselation (or triangulation in 2D), which was introduced by Voronoi [20] and extended by Delaunay [21]. The Delaunay triangles are constructed connecting the nodes whose Voronoi cells have common boundaries. From now on, we use 2D terminology although the discussion is also valid in 3D. The natural neighbours of any node are those defining its neighbouring Voronoi cells, that is, those connected with the node under consideration by the side of a Delaunay triangle.

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

## 2.2. Natural neighbour shape functions

In order to quantify the neighbour relation for any point $\mathbf{x}$ introduced into the tessellection, Sibson [15] used the concept of second-order Voronoi cells $T_{ij}$. $T_{ij}$ is defined as the locus of the points $\mathbf{x}$ such that its closest node is $n_i$ and the second closest node is $n_j$, i.e.

$$T_{ij} = \{ \mathbf{x} \in \mathbb{R}^n : d(\mathbf{x}, \mathbf{x}_i) < d(\mathbf{x}, \mathbf{x}_j) < d(\mathbf{x}, \mathbf{x}_k) \quad \forall k \neq i, \forall k \neq j, i \neq j \} \quad (2)$$

Let $\kappa$ be a Lebesgue measure (length, area or volume in 1D, 2D or 3D, respectively). Let $\kappa_x$ be the Lebesgue measure of $T_x$ and $\kappa_{xi}$ the Lebesgue measure of $T_{xi}$. The natural neighbour co-ordinates of $\mathbf{x}$ with respect to a natural neighbour node $n_i$ is defined as the ratio of the Lebesgue measure of second-order cell $T_{xi}$ to the Lebesgue measure of first-order cell $T_x$:

$$\phi_i(\mathbf{x}) = \frac{\kappa_{xi}}{\kappa_x} \quad (3)$$
<figure><img src="image_3.png" /><figcaption>Figure 2. Construction of natural neighbour shape functions.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Support of the natural shape function related to node $n_i$.</figcaption></figure>

If we consider the two-dimensional case shown in Figure 2, the natural neighbour co-ordinate or natural neighbour shape function is given by

$$ \phi_i(\mathbf{x}) = \frac{\text{Area}(\text{afghe})}{\text{Area}(\text{abcde})} \qquad (4) $$

Braun and Sambridge [12] have extended the Sibson theory [15] to deduce expressions for the derivatives of the interpolation functions, which can be applied in a Galerkin framework for solving partial differential equations, where trial and test functions are interpolated using natural neighbour shape functions. These authors refer to this method as the NEM.

Thus, the interpolation of a vector-valued function $\mathbf{u}^h(\mathbf{x})$ results

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

where $\mathbf{u}_i$ ($i \in [1, \dots, n]$) represents $\mathbf{u}(\mathbf{x}_i)$ ($\mathbf{x}_i$ are the co-ordinates of node $n_i$), $n$ the number of natural neighbour nodes of the point $\mathbf{x}$, and $\phi_i(\mathbf{x})$ is the shape function related to node $n_i$ and given by Equation (4).

It turns out that the support of $\phi_i(\mathbf{x})$ is the union of the $n$ circles (spheres in 3D) passing through the vertices of the $n$ Delaunay triangles (tetrahedra) containing the node $n_i$ (in this case $n$ is the number of natural neighbours of node $n_i$). The support of a node $n_i$ in a particular nodal distribution is depicted in Figure 3.
## 2.3. Natural element interpolation properties

In this section, the main properties of the natural neighbour shape functions are briefly reviewed. See Reference [13] for more details.

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

$$ \phi_i(\mathbf{x}_j) = \delta_{ij} \qquad (6) $$

which implies that $\mathbf{u}^h(\mathbf{x}_i) = \mathbf{u}_i$. Thus, the imposition of the essential boundary conditions is direct for the concerned nodes, in contrast to most of the other meshless approximations where $\mathbf{u}^h(\mathbf{x}_i) \neq \mathbf{u}_i$. This property and the linear variation of the NEM approximation on convex boundaries allows to easily impose essential boundary conditions.

* **The partition of unity:** By construction we have the following relation:

$$ \sum_{i=1}^{n} \phi_i(\mathbf{x}) = 1, \quad \forall \mathbf{x} \in \Omega \qquad (7) $$

This property is equivalent to a zeroth-order consistency, which means that the interpolant can exactly reproduce constant functions.

* **Linear consistency:** The natural neighbour shape functions also satisfy the local co-ordinate property [15], namely,

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

which combined with Equation (7), implies that the natural neighbour interpolant spans the space of linear polynomials (linear completeness).

* **Linear variation on convex boundaries:** The trial function is strictly linear on a Delaunay triangle edge located on the boundary of the convex hull. The proof of this property can be found in Reference [13].

* **Smoothness:** Natural neighbour shape functions are $C^\infty$ at any point except at the nodes, where they are only $C^0$, and on the boundary of the Delaunay circles (Figure 3) (spheres in 3D) where they are only $C^1$, because of the discontinuity on the number of natural neighbours nodes across these boundaries. Hiyoshi and Sughara [22] have shown that the Sibson interpolant belong to a more general class of Voronoi-based interpolants, called kth-order standard co-ordinates, proving that the interpolant generated by the kth-order standard co-ordinates have $C^k$ continuity on the Delaunay circles (spheres) boundaries. In this context, the Sibsonian and non-Sibsonian (Laplace) co-ordinates [23] results to be the standard co-ordinates of order 1 and 0, respectively.

## 2.4. Issues on non-convex domains

Since the natural neighbour interpolant is constructed on the basis of the Voronoi diagram and Delaunay tessellation, issues occur when non-convex domains are considered. The Voronoi diagram is constructed filling the convex hull of the cloud of nodes, without taking into account the actual domain boundary.

Unfortunately, on non-convex domains, some nodes lying on non-convex boundaries can share a Voronoi cell facet and thus they become natural neighbours. In this situation, the properties
<figure><img src="image_5.png" /><figcaption>Figure 4. Influences from nodes lying on non-convex boundaries.</figcaption></figure>

cited above are not satisfied anymore. The lack of linear precision on non-convex boundaries is also a consequence of this problem. To illustrate this fact, we consider the non-convex domains depicted in Figure 4, where a point **x** is introduced. In order to compute the interpolation at that point a new Voronoi cell is created, whereas the Voronoi cells related to the neighbour nodes of point **x** are modified (see Figure 4). We can see in Figure 4 (left) that the Voronoi cell just created shares edges with the Voronoi cells associated to nodes $n_1, n_2, n_3$ and $n_4$ that then result natural neighbours of point **x**. However, in order to impose the essential boundary condition at point **x** we require both, the shape function associated with node $n_3$ must vanish at point **x** and the influence of node $n_4$ on the solution at point **x** must be removed, i.e. $n_4$ cannot be a natural neighbour of point **x**.

Several techniques have been proposed to overcome this problem lying to the non-convex domains. A coupling between the finite element and natural element method [23] or the use of the $\alpha$-shape concept [24] can be cited among many other possibilities.

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

### 3. THE CONSTRAINED NATURAL ELEMENT METHOD (C-NEM)

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

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

Organ et al. have introduced in Reference [25] a visibility criterion to define the domain of influence of a node near a crack tip and other non-convex domains. In their approach, the point $\mathbf{x}$ is in the domain of influence of a node $n_i$ if $\mathbf{x}$ is within the region where the shape function associated with $n_i$ is non-zero ($\phi_i > 0$) and it is visible from node $n_i$ when the boundaries are assumed opaque. Thus, we can define

$$ \Omega_i = \{\mathbf{x} | \mathbf{x} \in \Omega, \phi_i(\mathbf{x}) > 0\} \qquad (9) $$

$$ \Omega_i^{VIS} = \{\mathbf{x} | S_{\mathbf{x}} \to \mathbf{x}_i \cap \Gamma = \emptyset\} $$

where $S_{\mathbf{x} \to \mathbf{x}_i}$ represents the straight line connecting $\mathbf{x}$ and $n_i$, i.e.

$$ S_{\mathbf{x} \to \mathbf{x}_i} = \mathbf{x} + \lambda(\mathbf{x}_i - \mathbf{x}), \quad \lambda \in [0, \dots, 1] \qquad (10) $$

where $\mathbf{x}_i$ is the vector defining the position of node $n_i$ and $\Gamma$ is the domain boundary (including interior parts). Thus, the domain of influence of a node $n_i$, $\Omega^*$, is defined by

$$ \Omega_i^* = \Omega_i^{VIS} \cap \Omega_i \qquad (11) $$

### 3.2. C-NEM interpolation

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

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

where $V$ is the number of natural neighbours visible from point $\mathbf{x}$ and $\phi_i^C$ is the constrained natural neighbour shape function. The computation of the C-n-n shape function is similar to natural neighbour shape function, when one proceeds on the constrained Voronoi diagram introduced later.

### 3.3. Support of constrained natural neighbour shape functions

In the C-NEM, the support of $\phi_i^C$ consists of the union of the $n$ circles (spheres in 3D), passing through the vertices of the Delaunay triangles (tetrahedra in 3D) defined by the node under consideration $n_i$, and the $V$ natural neighbours visible from node $n_i$. Furthermore, as stated by Edelsbrunner [26], these circles can extend across the boundaries, but cannot contain non-visible neighbours. In Figure 5, the intersection between the support of a constrained natural neighbour shape function $\phi_i^C$ and the boundary is depicted.

The discontinuity is explicitly described by nodes. Owing to the salient shape of the support, which passes through the visible neighbour nodes, the shape functions are strictly continuous near a crack tip, without requiring the use of particular techniques to ensure the continuity [25].
<figure><img src="image_6.png" /><figcaption>Figure 5. Supports of the natural neighbour and constrained natural neighbour shape functions.</figcaption></figure>

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

### 3.4. Constrained Voronoi diagram

In this section, the constrained Voronoi diagram (CVD) is introduced to provide a suitable background for the computation of natural neighbour shape function over any kind of geometry (convex or not). The constrained Voronoi diagram (or extended) Voronoi diagram, is introduced by Seidel in Reference [27]. In our work, we prefer the term constrained instead of extended to avoid confusions with the extended finite element terminology [4]. We consider the definition of the CVD provided in Reference [26]: Imagine $\mathbb{R}^2$ is a sheet of paper $\Sigma_0$, with the points of the set of nodes $S$ and the line segments defining the boundary in a set $L$, drawn on it. For each $l_i \in L$, we cut $\Sigma_0$ open along $l_i$ and glue another sheet $\Sigma_i$, which also cut open along $l_i$. The gluing is done around $l_i$ such every traveller who cross $l_i$ switches from $\Sigma_0$ to $\Sigma_i$ and vice versa. A schematic view of the particular gluing necessary to achieve that effect is illustrated in Figure 6(a). We know what it means for two points on the primary sheet to be visible from each other. For other pairs we need a more general definition. For $i \neq 0$, points $x_0 \in \Sigma_0$ and $y_i \in \Sigma_i$ are visible if $xy$ crosses $l_i$, and $l_i$ is the first constraining line segment crossed if we traverse $xy$ in the direction from $x$ to $y$.

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

$$T_i^C = \{ \mathbf{x} \in \mathbb{R}^n: d(\mathbf{x}, \mathbf{x}_i) < d(\mathbf{x}, \mathbf{x}_j), \forall j \neq i, S_{\mathbf{x} \to n_i} \cap \Gamma = \emptyset, S_{\mathbf{x} \to n_j} \cap \Gamma = \emptyset \} \quad (13)$$

where $\Gamma$ is the domain boundary, composed with segments $l_i \in L$.

A generalization of the constrained Delaunay triangulation to 3D does not exist without adding new nodes, as shown in Reference [28]. Nevertheless, some techniques to construct 3D constrained Delaunay tesselations are provided in References [29, 30] by addition of Steiner points.
<figure><img src="image_7.png" /><figcaption>Figure 6. Interpretation of the constrained Voronoi diagram: (a) schematic view of the CVD; and (b) intersection between CVD and the domain closure.</figcaption></figure>

The constrained Voronoi diagram has been introduced for two reasons:

(i) Once such diagram is constructed, classical algorithms for the computation of the shape functions, Watson-Bowyer [12] or Lasserre [31] can be applied directly because connections between non-visible natural neighbours in non-convex domains have been removed.

(ii) Recently a new integration method, the SCNI [18] has been introduced in the meshless methods framework for improving both the accuracy and the efficiency of the numerical integration of the weak formulation. This technique, as proved later, is specially appropriate in the context of NEM, because the integration cells match the constrained Voronoi cells.

### 3.5. Computation of the C-n-n shape functions

Once the constrained Voronoi diagram is computed, selection of constrained natural neighbours is a direct task because this diagram only contains connections between mutually visible nodes.

In Figure 7(a), an example of the computation of a natural shape function making use of the standard Voronoi diagram is presented. We can notice that the $n_1, n_2, n_3$ and $n_4$ nodal values will contribute to the interpolation in **x** (located on the non-convex domain boundary). We can observe that the computed solution in **x** is not good because its value does not result from a linear interpolation between the nodal values $n_1$ and $n_3$. Thus, the imposition of an essential boundary condition is disturbed by the influence of internal ($n_2$) and external and consequently non-visible ($n_4$) nodes.

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

We can summarize these facts in the following manner: the second-order Voronoi cells related to nodes located on the domain boundary (convex or not), have an infinite area because
<figure><img src="image_8.png" /><figcaption>Figure 7. Computation of natural neighbour shape functions over a non-convex domain: (a) over standard Voronoi diagram; and (b) over constrained Voronoi diagram.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. View showing the linearity of the C-NEM shape function associated with node $n_i$ along a discontinuity or a non-convex boundary.</figcaption></figure>

the influence between non-visible nodes has been removed during the construction of the constrained Voronoi diagram. In that case, influences of interior nodes also vanish. Thus, essential boundary conditions can be easily enforced on any domain boundary (convex or not), where the approximation is in fact linear (Figure 8).
### Remark 1

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

### Remark 2

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

# 4. NUMERICAL EXAMPLES

## 4.1. Elastostatic formulation

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

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

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

The constitutive relation is given by

$$ \sigma = C\epsilon \qquad (15) $$

where $\sigma$ and $\epsilon$ are the vector form of the stress and linearized strain tensor (symmetric part of the gradient of displacements tensor), respectively, which are written using the conventional notations:

$$ \sigma = \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{bmatrix}, \quad \epsilon = \begin{bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ 2\epsilon_{12} \end{bmatrix} \qquad (16) $$

and $C$ is the elastic tensor matrix form.

The essential and natural boundary conditions are

$$ \mathbf{u} = \bar{\mathbf{u}} \quad \text{on } \Gamma_u \qquad (17) $$

$$ \sigma \mathbf{n} = \bar{\mathbf{t}} \quad \text{on } \Gamma_t $$

where $\Gamma = \Gamma_u \cup \Gamma_t$ is the boundary of $\Omega$, $\mathbf{n}$ is the unit outward vector defined on $\Gamma$, and $\bar{\mathbf{u}}$ and $\bar{\mathbf{t}}$ the prescribed displacements and tractions, respectively.
The variational formulation (principle of virtual work) 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} \epsilon^* \cdot \sigma \, d\Omega = \int_{\Omega} v^* \cdot \mathbf{b} \, d\Omega + \int_{\Gamma_t} v^* \cdot \bar{\mathbf{t}} \, d\Gamma, \quad \forall v^* \in H_1^0(\Omega) \quad (18) $$

where $H_1(\Omega)$ and $H_1^0(\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 $v^*$, the following linear system of equations is obtained after numerical integration:

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

where $\mathbf{d}$ is the vector containing the nodal displacements, and the matrix $\mathbf{K}$ is given by

$$ \mathbf{K} = \int_{\Omega} \mathbf{B}' \mathbf{C} \mathbf{B} \, d\Omega \quad (20) $$

where $\mathbf{C}$ is the usual matrix form of the elastic tensor

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

$\mathbf{N}$ is the matrix containing the shape functions:

$$ \mathbf{N} = \begin{bmatrix} \phi_1 & 0 & \phi_2 & 0 & \dots & \phi_N & 0 \\ 0 & \phi_1 & 0 & \phi_2 & \dots & 0 & \phi_N \end{bmatrix} \quad (22) $$

$\mathbf{B}$ is the matrix containing the shape functions derivatives:

$$ \mathbf{B} = \begin{bmatrix} \phi_{1,x} & 0 & \phi_{2,x} & 0 & \dots & \phi_{N,x} & 0 \\ 0 & \phi_{1,y} & 0 & \phi_{2,y} & \dots & 0 & \phi_{N,y} \\ \phi_{1,y} & \phi_{1,x} & \phi_{2,y} & \phi_{2,x} & \dots & \phi_{N,y} & \phi_{N,x} \end{bmatrix} \quad (23) $$

The constitutive matrix $\mathbf{C}$ in isotropic linear elasticity results for plane strain:

$$ \mathbf{C} = \frac{E}{(1 - 2\nu)(1 + \nu)} \begin{bmatrix} 1 - \nu & \nu & 0 \\ \nu & 1 - \nu & 0 \\ 0 & 0 & \frac{1 - 2\nu}{2} \end{bmatrix} \quad (24) $$
and in the case of plane stress:

$$ C = \frac{E}{(1 - \nu^2)} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} \qquad (25) $$

## 4.2. Stabilized numerical integration

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

In this work, we use the integration procedure proposed by Chen et al. [18]: the stabilized conforming nodal integration (SCNI). In this approach, a strain smoothing stabilization is performed to stabilize the nodal integration.

The SCNI procedure is based on the assumed strain method, in which a modified gradient is introduced at the integration point (node) [18]:

$$ \tilde{\nabla} u_j^h(\mathbf{x}_k) = \frac{1}{A_k} \int_{\Omega_k} \nabla u_j^h(\mathbf{x}_k) \, d\Omega \qquad (26) $$

where $\mathbf{x}_k$ are the co-ordinates of node $n_k$.

Thus, the modified strain vector is given by

$$ \tilde{\epsilon}^h(\mathbf{x}_k) = \frac{1}{A_k} \int_{\Omega_k} \epsilon^h(\mathbf{x}) \, d\Omega = \frac{1}{A_k} \int_{\Omega_k} \begin{bmatrix} \frac{\partial u^h}{\partial x} \\ \frac{\partial v^h}{\partial y} \\ \frac{\partial u^h}{\partial y} + \frac{\partial v^h}{\partial x} \end{bmatrix} \, d\Omega \qquad (27) $$

Now, by applying the divergence theorem, the following equation is deduced:

$$ \tilde{\epsilon}^h(\mathbf{x}_k) = \frac{1}{A_k} \int_{\Gamma_k} \begin{bmatrix} u^h n_x \\ v^h n_y \\ u^h n_y + v^h n_x \end{bmatrix} d\Gamma \qquad (28) $$

where $\Omega_k$ is a representative domain related to node $n_k$ and $\Gamma_k$ and $A_k$ are its boundary and area, respectively, as depicted in Figure 9.

Introducing now the natural element approximation, i.e. $u^h$ and $v^h$ according to Equation (12) we obtain

$$ \tilde{\epsilon}^h(\mathbf{x}_k) = \tilde{\mathbf{B}}_k \mathbf{d} \qquad (29) $$
<figure><img src="image_10.png" /><figcaption>Figure 9. Integration domains for the SCNI procedure: (a) intersection between the Voronoi diagram and the domain closure; and (b) constrained and bounded Voronoi diagram.</figcaption></figure>

which is expressed explicitly by

$$ \tilde{c}^h(\mathbf{x}_k) = \begin{bmatrix} \tilde{b}_{x1}(\mathbf{x}_k) & 0 & \tilde{b}_{x2}(\mathbf{x}_k) & 0 & \dots & \tilde{b}_{xN}(\mathbf{x}_k) & 0 \\ 0 & \tilde{b}_{y1}(\mathbf{x}_k) & 0 & \tilde{b}_{y2}(\mathbf{x}_k) & \dots & 0 & \tilde{b}_{yN}(\mathbf{x}_k) \\ \tilde{b}_{y1}(\mathbf{x}_k) & \tilde{b}_{x1}(\mathbf{x}_k) & \tilde{b}_{y2}(\mathbf{x}_k) & \tilde{b}_{x2}(\mathbf{x}_k) & \dots & \tilde{b}_{yN}(\mathbf{x}_k) & \tilde{b}_{xN}(\mathbf{x}_k) \end{bmatrix} \begin{bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ \vdots \\ u_N \\ v_N \end{bmatrix} = \tilde{\mathbf{B}}_k \mathbf{d} \qquad (30) $$

where $\tilde{b}_{xj}(\mathbf{x}_k)$ and $\tilde{b}_{yj}(\mathbf{x}_k)$ are defined by

$$ \tilde{b}_{xj}(\mathbf{x}_k) = \frac{1}{A_k} \int_{\Gamma_k} \phi_j(\mathbf{x}) n_x(\mathbf{x}) \, d\Gamma \qquad (31) $$

$$ \tilde{b}_{yj}(\mathbf{x}_k) = \frac{1}{A_k} \int_{\Gamma_k} \phi_j(\mathbf{x}) n_y(\mathbf{x}) \, d\Gamma \qquad (32) $$

Obviously, a lot of components in the matrix $\tilde{\mathbf{B}}_k$ are zero due to the compact support of the shape functions. Chen et al. [18] propose to use the intersection between the Voronoi
<figure><img src="image_11.png" /><figcaption>Figure 10. (a) Model geometry of bimaterial cylinder; and (b) constrained Voronoi diagram of bimaterial cylinder with irregular spaced nodes.</figcaption></figure>

diagram and the domain closure to define the integration cells $\Omega_k$. Introducing the C-NEM approximation, the smoothed strain and the nodal integration, the global stiffness matrix is obtained by assembling the contribution of each node $n_k$:

$$ \mathbf{K} = \sum_k \tilde{\mathbf{K}}_k = \sum_k A_k \tilde{\mathbf{B}}_k^t \mathbf{C} \tilde{\mathbf{B}}_k \quad (33) $$

The stress in each cell is then deduced from

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

The application of the divergence theorem avoid the computation of the shape functions derivatives, and only the evaluation of the shape functions on the boundary of the Voronoi cells is required. The stabilized nodal integration has been successfully applied by Gonzalez et al. [19] in the context of the $\alpha$-NEM, where a significant improvement in the accuracy was noticed.

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

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

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

$$u_r = ar + \frac{b}{r} \qquad (35)$$

In our case we can write the displacements in each domain in the form

$$u_r = a^1 r + \frac{b^1}{r} \quad R_{\text{int}} \le r \le R_1 \qquad (36)$$

and

$$u_r = a^2 r + \frac{b^2}{r} \quad R_1 \le r \le R_{\text{ext}} \qquad (37)$$

from which the strains and stresses can be deduced

$$\begin{align}
\epsilon_{rr}^1 &= a^1 - \frac{b^1}{r^2} & \epsilon_{rr}^2 &= a^2 - \frac{b^2}{r^2} \nonumber \\
\epsilon_{00}^1 &= a^1 + \frac{b^1}{r^2} & \epsilon_{00}^2 &= a^2 + \frac{b^2}{r^2} \tag{38} \\
\epsilon_{r0}^1 &= 0 & \epsilon_{r0}^2 &= 0 \nonumber
\end{align}$$

$$\begin{align}
\sigma_{rr}^1 &= a^1(\lambda^1 + 2\mu^1) - \frac{2b^1\mu^1}{r^2} & \sigma_{rr}^2 &= a^2(\lambda^2 + 2\mu^2) - \frac{2b^2\mu^2}{r^2} \nonumber \\
\sigma_{00}^1 &= a^1(\lambda^1 + 2\mu^1) + \frac{2b^1\mu^1}{r^2} & \sigma_{00}^2 &= a^2(\lambda^2 + 2\mu^2) + \frac{2b^2\mu^2}{r^2} \tag{39} \\
\sigma_{r0}^1 &= 0 & \sigma_{r0}^2 &= 0 \nonumber
\end{align}$$

where $\lambda^i$ and $\mu^i$, $i = 1, 2$, denote the Lame's coefficients of both materials, related to elasticity modulus and Poisson's coefficient by

$$\mu = \frac{E}{2(1 + \nu)} \quad \text{and} \quad \lambda = \frac{\nu E}{(1 - 2\nu)(1 + \nu)} \qquad (40)$$

The parameters $a^i$ and $b^i$, $i = 1, 2$, can be computed by applying the following boundary and transmission conditions assuming plane stress:

$$\begin{align}
-\sigma_{rr}(R_{\text{int}}) &= P \nonumber \\
\sigma_{rr}(R_{\text{ext}}) &= 0 \nonumber \\
\sigma_{rr}^{1}(R_{1}) &= \sigma_{rr}^{2}(R_{1}) \nonumber \\
u_{r}^{1}(R_{1}) &= u_{r}^{2}(R_{1}) \tag{41}
\end{align}$$

where P represents the internal pressure.

Numerical results of displacements and stress are compared to the analytical solution for a cloud of 315 regularly distributed nodes in Figures 11 and 12 ($R_{\text{int}} = 0.5$ mm, $R_1 = 0.7$ mm, $R_{\text{ext}} = 0.9$ mm, $v_1 = 0.3$, $E_1 = 210000$ MPa, $v_2 = 0.2$, $E_2 = 105000$ MPa, $P = 1$ MPa). The C-NEM is in good agreement with the analytical solution, and the $\sigma_{00}$ discontinuity across the
<figure><img src="image_12.png" /><figcaption>Figure 11. Displacements across the interface of the bimaterial cylinder.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12. Stresses across the interface of the bimaterial cylinder.</figcaption></figure>

interface is accurately reproduced. Figure 13 and Table I group the errors using the energy norm for the different clouds of nodes and both regular and irregular nodal distributions, where $\tilde{\sigma}^h$ is given by Equation (34). In Figure 13, $h$ denotes the radial distance between two nodes.

### 4.4. A second example: a mode-I crack analysis

A crack analysis using the NEM can be found in Reference [13], where due to the symmetry of the problem the computation domain is actually convex. It shows that the NEM solution convergence is faster than the finite element one. In the present example, a mode-I cracked plate is considered. The geometry of the model is depicted in Figure 14(a).
<figure><img src="image_14.png" /><figcaption>Figure 13. Rate of convergence of the C-NEM for the bimaterial cylinder problem using a regular nodal distribution.</figcaption></figure>

<table><caption>Table I. Error in the energy norm for the bicylinder problem.</caption><thead><tr><th rowspan="2">Number of nodes</th><th>||σ̃<sup>h</sup> - σ<sub>ex</sub>||<sub>E</sub></th><th>||σ̃<sup>h</sup> - σ<sub>ex</sub>||<sub>E</sub></th></tr><tr><th>Regular nodal distribution</th><th>Irregular nodal distribution</th></tr></thead><tbody><tr><td>147</td><td>1.13 × 10<sup>-4</sup></td><td>1.21 × 10<sup>-4</sup></td></tr><tr><td>315</td><td>7.41 × 10<sup>-5</sup></td><td>8.10 × 10<sup>-5</sup></td></tr><tr><td>687</td><td>4.95 × 10<sup>-5</sup></td><td>5.48 × 10<sup>-5</sup></td></tr><tr><td>1863</td><td>3.06 × 10<sup>-5</sup></td><td>3.46 × 10<sup>-5</sup></td></tr></tbody></table>

The analytical solution in the neighbourhood of the crack tip is given by [33]

$$ \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 (42) $$

$$ \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 (43) $$

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

In order to investigate the convergence of the numerical solution in the whole domain we enforce the boundary conditions $F_d = \sigma_{ex}n$ on the boundary of the problem domain.

In this example, $L = 1$ mm, $a = L/2$ mm and $K_I = 1$. In order to illustrate the ability of the C-NEM to handle the crack behaviour, no symmetry was considered.
<figure><img src="image_15.png" /><figcaption>Figure 14. (a) Geometry for the mode-I crack problem; and (b) exact boundary conditions.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. Convergence of the C-NEM solution for the mode-I crack problem.</figcaption></figure>

The convergence of the C-NEM solution is depicted in Figure 15. The parameter *h* denotes the distance between two nodes in a regular set of nodes. We can notice that the convergence rate is 0.42, which is approximatively the expected convergence rate (0.5) obtained with triangular linear finite elements.

### 4.5. A third example: a plate containing an internal crack

In this example, we consider the problem defined in Figure 16 where $2a = 0.2$ mm, $2W = 2$ mm and $\sigma_{22}^{\infty} = 1$ MPa. Plane stress conditions are assumed with an elastic modulus of $E = 1$ MPa and a Poisson coefficient of $\nu = 0.3$. Firstly, we consider the case of $\gamma = 0$ whose exact solution in an infinite plate was obtained by Muskelishvili [34]. For $\gamma = 0$ and $\theta = 0$ in
<figure><img src="image_17.png" /><figcaption>Figure 16. Square plate with an inclined crack: (a) geometry and loads; and (b) distribution of nodes.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17. Zoom of the constrained Voronoi diagram in the neighbourhood of the crack (detail).</figcaption></figure>

Figure 16, this solution results

$$ \sigma_{22}(\theta=0, \gamma=0, r) = \sigma_{22}^{\infty} \frac{a+r}{\sqrt{r(2a+r)}} \quad (45) \\ \sigma_{11}(\theta=0, \gamma=0, r) = \sigma_{22}^{\infty} - \sigma_{22} $$

When the crack size is much smaller than the domain size, Muskelishvili's solution can be taken as a reference solution (Figure 17).

A cloud of 760 nodes is considered in the C-NEM simulation (400 nodes are uniformly distributed, other 20 nodes are located on each crack side and the remaining 320 nodes are concentrated in the crack tip zones) as depicted in Figure 16 without assuming any symmetry condition. Numerical results are reported in Figure 18.
<figure><img src="image_19.png" /><figcaption>Figure 18. Radial stresses ahead of the crack tip for $\gamma=0$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19. Stress intensity factor $K_I$ for different values of the crack orientation $\gamma$.</figcaption></figure>

Now, we consider the case of an inclined crack. The stress intensity factor (SIF) $K_I$ has been computed using the contour form of the interaction integral method [35], the different values of interest being easily interpolated anywhere in the domain. The size of the square domain used in the computation is $0.6a \times 0.6a$ (see Figure 6). Numerical results concerning the intensity factor $K_I$ for different values of $\gamma$ are compared in Figure 19 with the exact ones given in Reference [34]:

$$K_I = \sigma_{22}^{\infty} \sqrt{\pi a} \cos^2 \gamma \qquad (46)$$

Accuracy can be improved by increasing the nodal density in the crack tip neighbourhood. Another possibility is to extrinsically enhance the trial function by adding terms reproducing the local $\sqrt{r}$ behaviour. As the shape functions $\phi_i^C(x)$ define a partition of unity, the enrichment of the C-NEM trial function in crack problems can be achieved in the partition of unity framework
[4]. Let $\Omega_0 \subset \Omega$ be a sub-domain in a region around the crack tip where the crack tip singularity dominates. Then the enriched function for two-dimensional crack problems could follow the form suggested by Sukumar in Reference [13]:

$$ \mathbf{u}^h(\mathbf{x}) = \sum_{i=1}^{n} \phi_i^c(\mathbf{x})\mathbf{u}_i + \sqrt{r} \sum_{j=1}^{m} \phi_j^c(\mathbf{x})\mathbf{a}_j \quad (47) $$

where $\mathbf{a}_j$ are additional nodal coefficients that are associated with the $m$ nodes in $\Omega_0$.

## 5. CONCLUSION

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

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

(ii) The ability to reproduce discontinuities across interfaces.

(iii) Continuity inside the domain, which is not perturbed by the introduction of the visibility criterion.

(iv) Linear precision on any boundary (convex or not), without parameter dependence.

(v) A significant improvement of both efficiency and accuracy in the numerical integration by using the SCNI procedure, as demonstrated in Reference [19].

(vi) Linear consistency and partition of unity property at point of the domain. The second property allows functional enrichment.

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

## REFERENCES

1. Lucy LB. A numerical approach to the testing of fusion process. *The Astronomic Journal* 1977; **88**:1013–1024.

2. Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. *Computational Mechanics* 1992; **10**:307–318.

3. Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. *International Journal for Numerical Methods in Engineering* 1994; **37**:229–256.

4. Melenk JM, Babuška I. The partition of unity finite element method: basic theory and applications. *Computer Methods in Applied Mechanics and Engineering* 1996; **139**:289–314.

5. Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. *International Journal for Numerical Methods in Fluids* 1995; **21**:1081–1106.

6. Duarte CA, Oden JT. An H-p adaptative method using clouds. *Computer Methods in Applied Mechanics and Engineering* 1996; **139**:237–262.

7. Chen JS, Pan C, Wu CT, Liu WK. Reproducing kernel particle methods for large deformation analysis of nonlinear structures. *Computer Methods in Applied Mechanics and Engineering* 1996; **139**:195–229.
8. Günther FC, Liu WK. Implementation of boundary conditions for meshless methods. *Computer Methods in Applied Mechanics and Engineering* 1997; **163**:205-230.

9. Kaljevic I, Saigal S. An improved element free Galerkin formulation. *International Journal for Numerical Methods in Engineering* 1997; **40**:2953-2974.

10. Zhu T, Altruni SN. A modified collocation method and a penalty formulation for enforcing the essential boundary conditions in the element free Galerkin methods. *International Journal for Numerical Methods in Engineering* 1998; **21**:211-222.

11. Chen JS, Wang HP. New boundary conditions treatments in meshfree computation of contact problems. *Computer Methods in Applied Mechanics and Engineering* 2000; **187**:441-468.

12. Sambridge M, Braun J, McQueen M. Geophysical parameterization and interpolation of irregular data using natural neighbours. *Geophysical Journal International* 1995; **122**:837-857.

13. Sukumar N, Moran B, Belytschko T. The natural elements method in solid mechanics. *International Journal for Numerical Methods in Engineering* 1998; **43**:839-887.

14. Cueto E, Cegoñino J, Calvo B, Doblaré M. On the imposition of essential boundary conditions in natural neighbour Galerkin methods. *Communications in Numerical Methods in Engineering* 2003; **19**:361-376.

15. Sibson R. A vector identity for the Dirichlet tesselations. *Mathematical Proceedings of the Cambridge Philosophical Society* 1980; **87**:151-155.

16. Yvonnet J, Ryckelynck D, Lorong P, Chinesta P. Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi contraint—Méthode des éléments C-Naturels. *Revue Européenne des éléments Finis* 2003; **12**(4):482-509.

17. Krysl P, Belytschko T. Element-free Galerkin: convergence of the continuous and discontinuous shape functions. *Computer Methods in Applied Mechanics and Engineering* 1997; **148**:257-277.

18. Chen JS, Wu CT, Yoon Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. *International Journal for Numerical Methods in Engineering* 2001; **50**:435-466.

19. Gonzalez D, Cueto E, Martinez MA, Doblaré M. Numerical integration in natural neighbour Galerkin methods. *International Journal for Numerical Methods in Engineering* 2004; accepted for publication.

20. Voronoi GM. Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Deuxième mémoire: recherches sur les parallélloèdres primitifs. *Journal fur die Reine und Angewandte Mathematik* 1908; **134**:198-287.

21. Delaunay B. Sur la sphère vide. A la mémoire de Georges Voronoi *Izvestiya Akademicheskikh i Esteticheskikh Nauk* 1934; **7**:793-800.

22. Hiyoshi H, Sugura K. Improving continuity of Voronoi-based interpolation over Delaunay spheres. *Computational Geometry* 2002; **22**:167-183.

23. Sukumar N, Moran B, Semenov Y, Belikov VV. Natural neighbour Galerkin methods. *International Journal for Numerical Methods in Engineering* 2001; **50**:207-219.

24. Cueto E, Doblaré M, Gracia L. Imposing essential boundary conditions in the natural elements method by means of density-scaled alpha-shapes'. *International Journal for Numerical Methods in Engineering* 2000; **49**:519-546.

25. Organ D, Fleming M, Terry T, Belytschko T. Continuous meshless approximations for nonconvex bodies by diffraction and transparency. *New Advances in Adaptive Computational Mechanics* 1996; **19**:1-11.

26. Edelsbrunner H. *Geometry and Topology for Mesh Generation*. Cambridge University Press: Cambridge, 2001; p30.

27. Seidel R. Constrained Delaunay triangulations and Voronoi diagrams with obstacles. In *1978-1988 Ten Years IIG*, 1988; 178-191.

28. Schönhardt E. Über die Zerlegung von dreieckspolyedern in tetraeder. *Mathematische Annalen* 1928; **98**: 309-312.

29. Shewchuck JR. Tetrahedral mesh generation by Delaunay refinement. *Proceedings of the 14th Annual Symposium on Computational Geometry*, Minneapolis, Minnesota, Association for Computing Machinery, June, 1998; 86-95.

30. Shewchuck JR. Sweep algorithms for constructing higher-dimensional constrained Delaunay triangulations. *Proceedings of the 16th Annual Symposium on Computational Geometry*, Hong-Kong, Association for Computing Machinery, June, 2000; 350-359.

31. Sambridge M, Braun J, McQueen M. Computational methods for natural neighbour interpolation in two and three dimensions. *Proceedings of the Seventh Biennial Conference on Computational Techniques and Applications (CTAC 95)*. Melbourne, Australia, R. L. May, A. K. Easton (eds). 1996; 685-692.
32. Dolbow J, Belytschko T. Numerical integration of the Galerkin weak forme in meshfree methods. *Computational Mechanics* 1999; **23**:219–230.

33. Perrone N, Liebowitz H, Mulville D, Plikey W. *Fracture Mechanics*. University Press of Virginia: Charlottesville, 1978.

34. Lemaitre J, Chaboche JL. *Mechanics of Solid Materials*. Cambridge University Press: Cambridge, UK, 1990.

35. Yau J-F, Wang S-S, Corten HT. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity. *Journal of Applied Mechanics* 1980; **47**:335–341.