# Enrichment and coupling of the finite element and meshless methods

Antonio Huerta, Sonia Fernández-Méndez

<figure><img src="image_1.png" /></figure>

## To cite this version:

Antonio Huerta, Sonia Fernández-Méndez. Enrichment and coupling of the finite element and meshless methods. International Journal for Numerical Methods in Engineering, 2000, 48 (11), pp.1615-1636. 10.1002/1097-0207(20000820)48:113.0.CO;2-S. http://dx.doi.org/10.1002/1097-0207(20000820)48:113.0.CO;2-S

HAL Id: hal-01428629

https://hal.science/hal-01428629

Submitted on 6 Jan 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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Enrichment and coupling of the finite element and meshless methods

Antonio Huerta and Sonia Fernández-Méndez

Departament de Matemàtica Aplicada III, E.T.S. de Ingenieros de Caminos, Canales y Puertos,
Universitat Politècnica de Catalunya, Campus Nord, E-08034 Barcelona, Spain

A mixed hierarchical approximation based on finite elements and meshless methods is presented. Two cases are considered. The first one couples regions where finite elements or meshless methods are used to interpolate: continuity and consistency is preserved. The second one enriches a finite element mesh with particles. Thus, there is no need to remesh in adaptive refinement processes. In both cases the same formulation is used, convergence is studied and examples are shown.

KEY WORDS: adaptivity; *h*-p refinement; finite element method; meshless method; mixed interpolation; convergence

## 1. INTRODUCTION

Meshless or particle methods such as *reproducing kernel particle methods* (RKPM) [1, 2] *element-free Galerkin* (EFG) [3–8] or *smooth particle hydrodynamics* (SPH) [9, 10] among others (see [11, 12] for a general presentation), have nowadays proven their applicability in computational mechanics. They do not require to generate a mesh (a connectivity matrix) and thus, they are specially suited for certain problems, for instance adaptive refinement computations or discontinuous field problems (i.e. crack propagation problems [6]). Moreover, the interpolation functions in meshless methods are particularly attractive in the presence of high gradients, concentrated forces, and large deformations. However, particle methods suffer from an important computational cost that reduces their range of practical (engineering) applicability.

On the other hand, from a practical point of view, finite elements are less costly, implement Dirichlet boundary conditions in a simple way (no need for Lagrange multipliers), and, above all, they are widely used and trusted by practitioners. However, the relative cost of the mesh generation
process is, for some problems, very large. In particular, the cost of remeshing in adaptive refinement problems is clearly not negligible.

Several authors have already proposed to use mixed finite elements and meshless interpolations. The objective is always to use the advantages of each method. Belytschko et al. [5] already show how to couple finite elements near the Dirichlet boundaries and element-free Galerkin in the interior of the computational domain. This simplifies considerably the prescription of essential boundary conditions. They do a mixed interpolation in the transition region: area where both finite elements and particles have an influence. This mixed interpolation requires the substitution of finite element nodes by particles and the definition of ramp functions. Thus, the transition is of the size of one finite element and the interpolation is linear. With the same objectives, Hegen [13] couples the finite element domain and the meshless region with Lagrange multipliers. Here a new formulation is proposed. It follows the ideas of Belytschko et al. [5], generalizes them for any order of interpolation, suppresses the ramp functions, and does not require the substitution of nodes by particles. That is, as many particles as needed can be added where they are needed independently of the adjacent finite element mesh. This is done preserving the continuity of the solution and enforcing uniform consistency for the mixed interpolation.

Liu et al. [14] proposed a mixed interpolation with other goals and different formulations. They suggest to enrich the finite element approximation with particle methods. In fact, the following adaptive process seems attractive: (1) compute an approximation with a coarse finite element mesh, (2) do an a posteriori error estimation and (3) improve the solution with particles without any remeshing process. Meshless methods are ideal for such a procedure.

In this paper we present a unified and general formulation for mixed interpolations in both cases (coupling and enrichment). The formulation is developed for the EFG method. However, its generalization to other particle methods is straightforward.

In the following sections the formulation is developed, the applicability conditions are discussed and the convergence properties are presented. Finally, several examples are presented to illustrate the advantages of such a mixed approximation.

## 2. FUNDAMENTALS OF MESHLESS APPROXIMATIONS

This section will not be devoted to develop or discuss meshless methods in detail or their relation with moving least-squares (MLS) interpolants. There are well-known references with excellent presentations of meshless methods, see for instance [11, 12, 15, 16]. Here some basic notions will be recalled in order to introduce the notation and the approach employed in the following sections.

Meshless methods, or particle methods, are based in a functional interpolation of the form

$$u(x) \approx u^{\rho}(x) = \sum_{j \in I^{\rho}} u(x_j) N_j^{\rho}(x) \quad (1)$$

given a number of particles {$x_j$}$_{j \in I^{\rho}}$ in the domain $\Omega$, $\Omega \subset \mathbb{R}^n$. The interpolation functions, $N_j^{\rho}(x)$, must be determined in a proper manner. In reproducing kernel particle methods [1, 2] (RKPM) the interpolation functions are obtained in the framework of the MLS interpolation. The element-free Galerkin method [3–8] (EFG) can be viewed as a particular case of the previous formulation [16].
Let us recall, in the context of the EFG method, how the interpolation (shape) functions are obtained. They are defined as

$$N_j^\rho(x) = \mathbf{P}^T(x_j) \alpha(x) \phi\left(\frac{x-x_j}{\rho}\right) \qquad (2)$$

where $\mathbf{P}^T(x) = \{p_0(x), p_2(x), \dots, p_l(x)\}$ includes a complete basis of the subspace of polynomials of degree $m$ and the vector $\alpha(x)$ in $\mathbb{R}^{l+1}$ is unknown. In one dimension, it is usual that $p_i(x)$ coincides with the monomials $x^i$, and, in this particular case, $l = m$.

The function $\phi(x)$ is a weighting function (positive, even and with compact support) which characterizes the meshless method. For instance, if $\phi(x)$ is continuous together with its first $k$ derivatives, the interpolation is also continuous together with its first $k$ derivatives. In RKPM, $\phi(x)$ is directly related to the window function of the reproducing kernel. In (2) the weighting function has been translated, centred in $x_j$, and its support scaled by the dilation parameter $\rho$.

The unknown vector $\alpha(x)$ is determined imposing the so-called reproducibility or consistency condition. It is, in fact, the MLS condition. This reproducibility condition imposes that $u^\rho(x)$ defined in (1) interpolates exactly polynomials of degree less or equal to $m$, i.e.

$$\mathbf{P}(x) = \sum_{j \in I^\rho} \mathbf{P}(x_j) N_j^\rho(x) \qquad (3)$$

For computational purposes, it is usual and preferable [16] to centre in $x_j$ and scale with $\rho$ also the polynomials involved in previous expressions. Thus, another expression for the shape functions is employed:

$$N_j^\rho(x) = \mathbf{P}^T\left(\frac{x-x_j}{\rho}\right) \alpha(x) \phi\left(\frac{x-x_j}{\rho}\right) \qquad (4)$$

which is similar to (2). The consistency condition becomes in this case

$$\mathbf{P}(0) = \sum_{j \in I^\rho} \mathbf{P}\left(\frac{x-x_j}{\rho}\right) N_j^\rho(x) \qquad (5)$$

which is equivalent to condition (3) when $\rho$ is constant everywhere. After substitution of (4) into (5) the usual linear system of equations, that determines $\alpha(x)$, is obtained

$$\mathbf{M}(x) \alpha(x) = \mathbf{P}(0) \qquad (6)$$

with

$$\mathbf{M}(x) = \sum_{j \in I^\rho} \mathbf{P}\left(\frac{x-x_j}{\rho}\right) \mathbf{P}^T\left(\frac{x-x_j}{\rho}\right) \phi\left(\frac{x-x_j}{\rho}\right) \qquad (7)$$

Notice that for each $x$ in $\mathbb{R}^n$ the previous sum only involves those particles whose support—the support of the weighting function $\phi$—includes $x$.

Section 3.2 presents a discussion on the necessary conditions for $\mathbf{M}(x)$ to be positive definite, namely, the requirements on the particle distribution and the value of the dilation parameter. This will allow to compute the vector $\alpha$ at each point and thus determine the shape functions, $N_j^\rho(x)$.
**Remark 1.** The consistency conditions (3) and (5) are equivalent if the dilation parameter $\rho$ is constant. When the dilation parameter varies at each particle another definition of the shape functions is recommended,

$$N_j^\rho(x) = \mathbf{P}^T \left( \frac{x-x_j}{\rho} \right) \alpha(x) \phi \left( \frac{x-x_j}{\rho_j} \right)$$

where $\rho_j$ is the dilation parameter associated to particle $x_j$, and a constant $\rho$ is employed in the scaling of the polynomials $\mathbf{P}$. Note that expression (4) is not directly generalized. The constant value $\rho$ is typically chosen as the mean value of all the $\rho_j$. The consistency condition in this case is also (5). It also imposes the reproducibility of the polynomials $\mathbf{P}$.

**Remark 2.** The dilation parameter $\rho$ characterizes the support of the shape functions $N_i^\rho(x)$. In fact, it plays a role similar to the element size in the finite element method. An $h$-refinement in finite elements can be produced in meshless methods decreasing the value of $\rho$ (this usually implies an increase in the number of particles). Liu et al. [16] proved convergence of the RKPM and, in particular, of EFG. The *a priori* bound is very similar to the bound in finite elements. The parameter $\rho$ plays the role of $h$, and $m$ (the order of consistency) plays the role of the degree of the interpolation polynomials in the finite element mesh.

**Remark 3.** Convergence properties depend on $m$ and $\rho$. They do not depend on the distance between particles because usually [10, 16] this distance is proportional to $\rho$, i.e. the ratio between the particle distance over the dilation parameter is of order one.

### 3. A HIERARCHICAL MIXED APPROXIMATION: FINITE ELEMENTS WITH EFG

Suppose, as discussed in the introduction, that the interpolation of $u(x)$ in $\Omega$, $\Omega \subset \mathbb{R}^n$, is done with both finite elements and EFG. The domain must include a set of nodes $\{x_i\}_{i \in I^h}$ with their associated shape functions $N_i^h(x)$, that are going to take care of the finite element contribution, $u^h(x)$, to $u(x)$, namely,

$$u^h(x) = \sum_{i \in I^h} u(x_i) N_i^h(x) \qquad (8)$$

There is also a set of particles $\{x_j\}_{j \in I^\rho}$ with their associated interpolation functions $N_j^\rho(x)$, that are going to take care of the meshless contribution,

$$u^\rho(x) = \sum_{j \in I^\rho} u(x_j) N_j^\rho(x) \qquad (9)$$

In the more general case, the domain $\Omega$ is the union of two non-disjoint subdomains, $\Omega = \Omega^h \cup \Omega^\rho$, where $\Omega^h$ denotes the subdomain where the $N_i^h$ have an influence in the approximation,

$$\Omega^h = \{x \in \Omega | \exists i \in I^h, N_i^h(x) \neq 0\}$$

and where $\Omega^\rho$ is the subdomain where at least one $N_j^\rho$ is non-zero,

$$\Omega^\rho = \{x \in \Omega | \exists j \in I^\rho, N_j^\rho(x) \neq 0\}$$
<figure><img src="image_3.png" /><figcaption>Figure 1. Coupled finite element and element-free Galerkin.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. Finite element enrichment with element-free Galerkin.</figcaption></figure>

In the region where only finite elements are present, $\Omega^h \setminus \Omega^\rho$, a standard, and thus consistent, finite element approximation is considered:

$$u(x) \approx u^h(x)$$

In the region where only particles have an influence, $\Omega^\rho \setminus \Omega^h$, the standard, and thus consistent, EFG approximation is considered:

$$u(x) \approx u^\rho(x)$$

However, in the area where both interpolations have an influence, $\tilde{\Omega} := \Omega^h \cap \Omega^\rho$, a mixed interpolation must be defined

$$u(x) \approx u^h(x) + u^\rho(x) \qquad (10)$$

The objective now is to develop a mixed functional interpolation, such as (10), with the desired consistency in $\tilde{\Omega}$, without any modification of the finite element shape functions $N_i^h$ and such that $u^\rho(x)$ is hierarchical. That is, an EFG contribution which should be zero at the finite element nodes, must be added to the standard finite element interpolation. Obviously, such a contribution must verify consistency conditions similar to those of standard meshless methods [2, 3, 16].

In the following sections this mixed interpolation is developed and discussed. In particular, the admissible particle distribution is detailed. Moreover, a priori convergence is studied when the number of particles is increased, when the number of nodes is increased, and when both particles and nodes are increased.

Moreover, two cases will be considered with the same formulation. The first one (coupled finite element and element-free Galerkin) requires that $\Omega^h \neq \Omega$ and $\Omega^\rho \neq \Omega$. That is, in a region of $\Omega$ only finite elements will be used, in another region only EFG are employed, and finally in a mixed area, $\tilde{\Omega}$, the solution is approximated using (10). In the second case, $\Omega = \Omega^h$ and $\tilde{\Omega} = \Omega^\rho$. That is, there is a complete finite element basis all over $\Omega$. Only in a reduced area, $\tilde{\Omega}$, particles are added to improve the interpolation (finite element enrichment with element-free Galerkin). Both situations are depicted in Figures 1 and 2, and developed in Sections 4 and 5.

### 3.1. Evaluation of the meshless shape functions $N_i^\rho$

In $\tilde{\Omega}$ the expression of the interpolation function is obtained after substitution of (8) and (9) into (10), namely

$$u(x) \approx \sum_{i \in I^h} u(x_i) N_i^h(x) + \sum_{j \in I^\rho} u(x_j) N_j^\rho(x) \qquad (11)$$
where $N^\rho(x)$ is defined, as previously, in (4), and, as before, the vector of unknown functions, $\alpha(x)$, is determined using the consistency condition. Now the reproducibility conditions impose that (11) must interpolate exactly a complete basis of polynomials of order less or equal to $m$. That is,

$$ \mathbf{P}(0) = \sum_{j \in I^\rho} \mathbf{P}\left(\frac{x-x_j}{\rho}\right) N_j^\rho(x) + \sum_{i \in I^h} \mathbf{P}\left(\frac{x-x_i}{\rho}\right) N_i^h(x) \quad (12) $$

which is the natural extension of (5). Note that Equation (12) can, if $\rho$ is constant everywhere, be rewritten as

$$ \mathbf{P}(x) = \sum_{j \in I^\rho} \mathbf{P}(x_j) N_j^\rho(x) + \sum_{i \in I^h} \mathbf{P}(x_i) N_i^h(x) \quad (13) $$

which shows more clearly the desired reproducibility condition. The linear system of equations that determines $\alpha$ is obtained once the definition of $N^\rho(x)$, Equation (4), is substituted in (12),

$$ \mathbf{M}(x) \alpha(x) = \mathbf{P}(0) - \sum_{i \in I^h} \mathbf{P}\left(\frac{x-x_i}{\rho}\right) N_i^h(x) \quad (14) $$

The least-squares matrix is identical to the matrix employed in the standard EFG method, Equation (7), namely

$$ \mathbf{M}(x) = \sum_{j \in I^\rho} \mathbf{P}\left(\frac{x-x_j}{\rho}\right) \mathbf{P}^T\left(\frac{x-x_j}{\rho}\right) \phi\left(\frac{x-x_j}{\rho}\right) $$

*Remark 4.* The particle shape functions $N_j^\rho$ are hierarchical. Note that at any node $x_k$, $k \in I^h$, the right-hand side of (14) is zero

$$ \mathbf{P}(0) - \sum_{i \in I^h} \mathbf{P}\left(\frac{x_k - x_i}{\rho}\right) N_i^h(x_k) = \mathbf{P}(0) - \sum_{i \in I^h} \mathbf{P}\left(\frac{x_k - x_i}{\rho}\right) \delta_{ik} = \mathbf{0} $$

Thus, the solution of (14) is trivial, $\alpha(x_k) = \mathbf{0}$. And therefore, from the definition of the particle shape functions, (4), it is easy to verify that the $N_j^\rho(x)$ are hierarchical, i.e. $N_j^\rho(x_k) = 0 \quad \forall j \in I^\rho, k \in I^h$.

## 3.2. Admissible particle distribution

As in standard EFG, matrix $\mathbf{M}(x)$ must be regular (invertible) everywhere, i.e. at each point $x \in \tilde{\Omega}$. Only the right-hand side of (14) differs from the EFG system of equations (6). Thus, as in EFG, the number of particles, their position and their related dilation parameters cannot be taken arbitrarily. In Liu et al. [16] there is an excellent definition of the admissible particle distribution. Here some essential details are recalled in order to discuss the validity of the mixed approximation.

At a point $x^* \in \tilde{\Omega}$, matrix $\mathbf{M}(x^*)$ can be viewed as a Gram matrix defined with the discrete scalar product

$$ \langle f, g \rangle_{x^*} = \sum_{j \in I^\rho} f\left(\frac{x^* - x_j}{\rho}\right) g\left(\frac{x^* - x_j}{\rho}\right) \phi\left(\frac{x^* - x_j}{\rho}\right) \quad (15) $$
<figure><img src="image_5.png" /><figcaption>Figure 3. Substitution of a finite element node by one particle. Non-admissible distribution.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4. Substitution of a finite element node by two particles. Admissible distribution.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Non-admissible distribution. $\tilde{\Omega}$ is under the influence of only one particle.</figcaption></figure>

and with the linear independent polynomials in $P(x)$. If the scalar product, $\langle \cdot, \cdot \rangle_{x^*}$, is degenerated the matrix $M(x^*)$ is singular. For instance, every point $x^* \in \tilde{\Omega}$ must lie in the area of influence of, at least, $l + 1 = \dim(M)$ particles. That is, the following condition is necessary:

$$ \mathrm{card} \left\{ x_j \mid j \in I^\rho, \phi\left(\frac{x^* - x_j}{\rho}\right) \neq 0 \right\} \geq \dim(M) = l + 1 $$

Moreover, the location of those particles is not arbitrary. In a $n$-dimensional space, i.e. $x^* \in \mathbb{R}^n$, the $n + 1$ particles needed for a linear interpolation must describe a non-degenerated $n$-simplex. For instance, in two dimensions ($n = 2$), $x^*$ must belong to the support of at least three shape functions associated to particles not aligned; or in three dimensions ($n = 3$), $x^*$ must belong to the support of at least four functions with particles not coplanar.

These restrictions are also valid for possible distributions of particles in a mixed interpolation. For instance, in a one-dimensional domain with an order one consistency (linear interpolation) a finite element node cannot be replaced by a single particle, see Figure 3. Two particles, with dilation parameters large enough, are needed in order to ensure that everywhere in $\tilde{\Omega}$ the scalar product, (15), does not degenerate. Figures 3 and 4 depict these situations. For each particle, its corresponding weighting function $\phi((x - x_i)/\rho)$ is plotted.

Figure 5 also shows a non admissible distribution of particles. In the region where both particle and finite element interpolations have an influence, $\tilde{\Omega}$, there are not enough particles (only one is present) to ensure the regularity of $M(x)$. An obvious solution for this problem, maintaining the same particle distribution, is to choose a dilation parameter large enough, see Figure 6.

**Remark 5.** The shape functions $N_j^\rho$ are hierarchical. Thus, in one dimension, the weighting functions $\phi((x - x_i)/\rho)$ can be truncated outside $\Omega^\rho$ and continuity of $N_j^\rho(x)$ is preserved, see
<figure><img src="image_8.png" /><figcaption>Figure 6. Approximation functions before and after imposing the consistency condition of order one.</figcaption></figure>

Figures 4 and 6. This property cannot be generalized to higher dimensions. In fact Sections 4 and 5 discuss this issue.

Remark 6. As previously indicated for EFG, the interpolation functions could be defined using (2) instead of (4). Thus, the scalar product needed to define $M(x^*)$ is now

$$ \langle f, g \rangle_{x^*} = \sum_{j \in I^\rho} f(x_j)g(x_j) \phi \left( \frac{x^* - x_j}{\rho} \right) \quad (16) $$

instead of (15). However, it is preferable to scale the polynomials $P(x)$ as done previously, see (4) and (15), because Gram matrices, such as $M(x^*)$, are easily ill-conditioned, specially with the trivial basis of polynomials. In general, with the translation to $x^*$ and the scaling with $\rho$, Gram matrices have lower condition numbers.

Remark 7. The aforementioned conditions for matrix $M(x)$ ensure its regularity. However, these conditions do not imply the solvability of the discrete variational problem. The ‘stiffness’ matrix for the global problem may be singular if a quadrature not accurate enough is employed. Moreover, as discussed in Section 5, Remark 9, when particles are added as an enrichment of finite elements, the shape functions associated to particles, $N^\rho$, are not linearly independent.

## 4. COUPLED FINITE ELEMENT AND ELEMENT-FREE GALERKIN

In this section, a new formulation, which generalizes the coupled formulation proposed by Belytschko et al. [5], is presented. This coupling between finite elements and EFG maintains both continuity and consistency everywhere, in particular, in the transition area. The major differences with the previously cited reference [5] are: (1) there is no need to replace nodes by particles, and (2) no ramp functions must be defined.

In fact, the generalization proposed here can be used for any order of consistency (it can go beyond linear elements and order one consistency). Moreover, this method allows to introduce as many particles as desired in the last element that defines the transition area, see Figures 7 and 8.

The computational domain $\Omega$ is divided into three non-disjoint regions: one where finite elements have an influence, $\Omega^h$, another where particles have an influence, $\Omega^\rho$, and finally, one region, $\tilde{\Omega}$, for the transition. In the latter, both particles and nodes define the interpolation, see Figure 1. Such a situation may be of interest if a computation with finite elements of degree $p$ needs to be refined in a region $\Omega^\rho$ without remeshing. The nodes of the original finite element mesh are removed in $\Omega^\rho$ but as many particles as needed are added in that region (see the crack propagation examples in the papers by Belytschko and co-workers [4–6, 8]).

It is important to notice that the approximation $u^h + u^\rho$ is continuous everywhere in $\Omega$ if the following conditions are met. First, the same order of consistency is imposed all over $\Omega$ (i.e. for both finite elements and particles), namely, $m = p$. And second, the domain of influence of
<figure><img src="image_9.png" /><figcaption>Figure 7. Coupled approximation functions with consistency of order one and two particles in the transition region $\tilde{\Omega}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Coupled approximation functions with consistency of order two and two different distribution of particles.</figcaption></figure>

particles, $\Omega^\rho$, coincides exactly with the region where finite elements do not have a complete basis. That is, no particles are added in 'complete' finite elements (i.e. elements where no node has been suppressed). Moreover, weighting functions $\phi$ are chopped off in those 'complete' finite elements, see Figure 6. In other words, $\Omega^\rho$ is the union of elements where at least one node has been removed.

The approximation $u^h + u^\rho$ is continuous as long as the shape functions $N_j^\rho$ are continuous. In spite of chopping off the weighting functions outside $\Omega^\rho$ the approximation maintains its regularity. This is due to the fact that $N_j^\rho(x) = 0$ over $\Omega^h \cap \partial\Omega^\rho$, with absolute independence of the fact that $\phi((x-x_j)/\rho) \neq 0$ over $\Omega^h \cap \partial\Omega^\rho$.

In $\Omega^h \setminus \Omega^\rho$ the finite element interpolation is complete and of order $m$. In particular, over $\Omega^h \cap \partial\Omega^\rho$ polynomials of degree less or equal to $m$ are interpolated exactly. Thus, it is easy to verify that

$$ \mathbf{P}(0) - \sum_{i \in I^h} \mathbf{P} \left( \frac{x - x_i}{\rho} \right) N_i^h(x) = \mathbf{0} \quad \text{over } \Omega^h \cap \partial\Omega^\rho $$

Recalling (14), the previous equation implies that $\alpha(x) = \mathbf{0}$, and consequently, $N_j^\rho(x) = 0$ for $x \in \Omega^h \cap \partial\Omega^\rho$, see Equation (4). Note that the previous rationale is independent of the spatial
<figure><img src="image_11.png" /><figcaption>Figure 9. Convergence of FEM and coupled FEM-EFG for a distribution of elements and particles shown in Figure 6.</figcaption></figure>

dimensions. Thus, in 1D, 2D or 3D, the shape functions $N_j^\rho$ are continuous independently of the truncation of the weighting functions $\phi((x - x_j)/\rho)$, see Figure 6.

Remark 8. In $\Omega^\rho \setminus \Omega^h$ finite elements have no influence,

$$ \mathbf{P}(0) - \sum_{i \in I^h} \mathbf{P} \left( \frac{x - x_i}{\rho} \right) N_i^h(x) = \mathbf{P}(0) $$

In this region, shape functions $N_j^\rho$ are identical to the standard EFG ones.

## 4.1. Convergence analysis

It is easy to verify that the mixed interpolation proposed preserves the convergence rate of FEM and EFG. Function

$$ u(x) = x^4 + 2x^3 \quad \text{is interpolated for } x \in \bar{\Omega} = [-1, 1] $$

The three regions of influence of finite elements, particles and the mixed interpolation are: $\Omega^h = [-h, 1]$, $\Omega^\rho = [-1, 0]$ and $\bar{\Omega} = [-h, 0]$, where $h$ is the size of finite elements, see Figure 6.

Figure 9 shows the convergence rate—logarithm of the error in $L^2([-1, 1])$ versus the total number of degrees of freedom—in two cases: standard linear finite elements and a coupled finite element-EFG approximation of order one. With this distribution of particles and with consistency of order one, this approach gives the same results as the one proposed by Belytschko et al. [5]. Similar conclusions can be drawn with other distributions of particles and order of consistency, see in the same figure the convergence results obtained with $m = 2$ and a particle distribution of Figure 8.
5. FINITE ELEMENT ENRICHMENT WITH ELEMENT-FREE GALERKIN

A finite element approximation can be improved (enriched) without any need of remeshing by adding particles. Particle methods have demonstrated their advantages in adaptive computations and their suitability to capture large gradients, concentrated loads and large deformations. Thus, enrichment of finite elements with meshless methods of the desired order seems an attractive option in these problems.

In this case, the region $\tilde{\Omega}$ where particles are added also maintains the original complete finite element interpolation, see Figure 2. In $\tilde{\Omega}$, the consistency of the mixed interpolation $m$ must be larger than the order of the finite element interpolation $p$. If consistency is set equal to $p$, finite elements can reproduce exactly polynomials up to degree $p$, thus

$$ \mathbf{P}(0) - \sum_{i \in I^h} \mathbf{P} \left( \frac{x - x_i}{\rho} \right) N_i^h(x) = \mathbf{0} \quad \forall x \in \tilde{\Omega} $$

and the solution of (14) is the trivial one, $\alpha = \mathbf{0}$. Consequently, the interpolation functions related to the particles $N_i^\rho$ are identically zero everywhere. Thus $\mathbf{P}(x)$ must include at least one polynomial not reproducible by the finite element interpolation, i.e. $m > p$.

As previously seen in Section 3 the shape functions $N_i^\rho$ are hierarchical. Thus, the interpolation is continuous in one-dimensional problems irrespective of the truncation of the weighting functions, $\phi(x)$ outside $\tilde{\Omega}$. In higher dimensions, continuity is not preserved as soon as the order of consistency is not constant and uniform everywhere in $\tilde{\Omega}$. In fact, the increase in consistency just mentioned in $\tilde{\Omega}$ will induce discontinuities in the approximation along $\partial\tilde{\Omega}$: functions $N_i^\rho$ are hierarchical but do not go to zero everywhere on $\partial\tilde{\Omega}$. If the approximation must be continuous a region surrounding $\tilde{\Omega}$ must be defined in which the interpolation functions $N_i^\rho$ go to zero with continuity. However, if $\partial\tilde{\Omega}$ coincides with an area where finite elements capture accurately the solution, those discontinuities due to the enrichment are going to be small.

*Remark 9.* Linear elements in 1D reproduce exactly polynomials of degree less or equal to one. In this case the first two equations of the system of Equations (14) are the consistency conditions:

$$ \sum_{j \in I^\rho} N_j^\rho(x) = 0 $$

$$ \sum_{j \in I^\rho} x_j N_j^\rho(x) = 0 $$

which correspond to the first two equations in (13). This implies that all the interpolation functions $N_j^\rho$ must verify these relations and, thus, they are no longer linearly independent. If every interpolation function is used in the resolution of the boundary value problem, the ‘stiffness’ matrix would be singular (two of its eigenvalues are zero). Thus, once the shape functions are evaluated, i.e. after (14) is solved, two of those interpolation functions are eliminated. Then, a linear set of interpolation functions is recovered and the ‘stiffness’ matrix remains regular. In general, it is necessary to suppress a $N_j^\rho$ (i.e. a particle) of the interpolation set for each polynomial in $\mathbf{P}(x)$ that finite elements are able to capture exactly.
<figure><img src="image_12.png" /><figcaption>Figure 10. Convergence for a mesh and meshless refinement: constant $h/\rho$ and $h \to 0$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Convergence for a mesh refinement: constant $\rho$ and $h \to 0$.</figcaption></figure>

## 5.1. Convergence analysis

A parametric analysis of convergence rates for this proposed method is shown next. The same function used previously is also used here,

$$u(x) = x^4 + 2x^3, \quad x \in \bar{\Omega} = [-1, 1]$$

with particles and finite elements everywhere. Finite elements are enriched everywhere adding particles and increasing the order of consistency. As before, $p$ is the degree of the finite element interpolation, and $m$ is the order of consistency obtained with the added particles. The increment of consistency $q$ is defined as

$$q := m - p$$

The error is evaluated in the $L^2(\Omega)$ norm. In Figure 10 the logarithm of the error is plotted against the logarithm of the number of degrees of freedom for different values of $p$ and $q$. Here, both finite element and meshless approximations are refined simultaneously (maintaining $h/\rho$ constant). Note that the order of the method is $O(h^{m+1})$. It is the same order that can be obtained with standard finite elements of degree $m = p+q$, or standard EFG with consistency of order $m$.

Figure 11 shows convergence results when the number of particles is kept constant but elements are refined. The order of the method is $O(h^{p+1})$ (identical to the order of finite elements alone)
<figure><img src="image_14.png" /><figcaption>Figure 12. Convergence for a meshless refinement: constant $h$ and $\rho \to 0$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13. Function $u(x)$ defined in (17).</figcaption></figure>

if the constant $\rho$ is large (four particles in $\Omega$). However, when the number of particles increases (256 particles in $\Omega$), i.e. small $\rho$, the order becomes $O(h^{m+1})$. Thus, if the density of particles is large, a refinement in finite elements of degree $p$ induces an order of convergence similar to finite elements of degree $m = p + q$.

Finally, Figure 12 shows the rate of convergence when refinement is only based on particles. That is, the finite element mesh is kept constant. If the element size, $h$, is small, the order observed is $O(\rho^q)$, but when the element size is large (four elements is $\Omega$) the mixed approximation does not converge as $\rho$ goes to zero.

*Remark 10.* Note that this convergence rates are also obtained with functions other than polynomials. In particular, similar results are obtained with the following function:

$$u(x) = \sin\left(\frac{7}{6}\pi(x+1)\right)\cos^3\left(\frac{35}{6}\pi(x+1)\right) \quad (17)$$

which is plotted in Figure 13.

In fact, this convergence analysis can further be exploited. An *a priori* error estimate can be obtained for the mixed approximation proposed in this section.
**Theorem 1.** Let $m$ be the order of consistency of the mixed approximation $u^h + u^\rho$, such that $m = p+q$, where $p$ is the order of the finite element approximation, $u^h$, and $q > 0$ is the order increment due to $u^\rho$. Suppose the following regularity conditions hold for the exact solution, $u$, and the weighting function, $\phi$: $u \in C^{m+1}(\bar{\Omega})$ and $\phi \in C^0(\bar{\Omega})$, where $\Omega$ is bounded and $\partial\Omega$ is smooth. Finally, assume that the element size is small enough, i.e.

$$ \frac{h}{\rho} \le \min_{p+1 \le r \le m} \left( \begin{array}{c} r \\ p+1 \end{array} \right)^{-1/(r-(p+1))} $$

Then,

$$ \|u - (u^h + u^\rho)\|_{L_\infty} \le h^{p+1} [C_1 h^q + C_2 \rho^q] |u|_{W_\infty^{m+1}} \quad (18) $$

where $C_1$ and $C_2$ are independent of the finite element size, $h$, and the dilation parameter, $\rho$, of the meshless approximation.

Note that on the left-hand side of the inequality the standard infinite norm over $\Omega$ is used, whereas on the right-hand side the seminorm $|\cdot|_{W_\infty^{m+1}}$ is employed. For the sake of clarity, multi-index notation is introduced: given the $n$-tuple $\alpha = (\alpha_1, \alpha_2, \dots, \alpha_n) \in \mathbb{N}^n$ and the non-negative integer $|\alpha| := \alpha_1 + \alpha_2 + \dots + \alpha_n$ then, by definition,

$$ |u|_{W_\infty^{m+1}} = \sum_{|\alpha|=m+1} \max_{x \in \bar{\Omega}} \left| \frac{\partial^{|\alpha|} u}{\partial x_1^{\alpha_1} \partial x_2^{\alpha_2} \dots \partial x_n^{\alpha_n}} \right| $$

It is important to remark that the error bound in (18) coincides with the convergence results shown in Figures 10–12. That is, when both $h$ and $\rho$ decrease simultaneously, the order of convergence is $p+q+1 = m+1$. When $h$ goes to zero while $\rho$ is kept constant, the order is either $p+1$ if $C_1 h^q < C_2 \rho^q$ or $m+1$ when $C_1 h^q \ge C_2 \rho^q$. And finally, convergence is ensured at a rate of $q$ when $\rho$ goes to zero provided that $C_1 h^q \le C_2 \rho^q$.

The previous theorem introduces a restriction on the element size which can be relaxed at a prize of obtaining a new error bound not as sharp.

**Theorem 2.** Under the same assumptions of Theorem 1 but with no restriction on the element size, the a priori error bound becomes

$$ \|u - (u^h + u^\rho)\|_{L_\infty} \le h^{p+1} [C_1 h^q + C_2 \rho^1] |u|_{W_\infty^{m+1}} $$

See Reference [17] for a detailed proof of the previous theorems. Moreover, following the ideas exposed in Liu et al. [16] Theorem 1 can be extended to the standard form in finite element analysis, see the proof in Reference [18].

**Theorem 3.** Let $m$ be the order of consistency of the mixed approximation $u^h + u^\rho$, such that $m = p+q$, where $p$ is the order of the finite element approximation, $u^h$, and $q > 0$ is the order increment due to $u^\rho$. Given $k$ such that $0 \le k \le p$, suppose the following regularity conditions hold for the exact solution, $u$, and the weighting function, $\phi$: $u \in C^{m+1}(\bar{\Omega})$ and $\phi \in C^k(\bar{\Omega})$, where
$\Omega$ is bounded and $\partial\Omega$ is smooth. Finally, assume that the element size is small enough, i.e.

$$ \frac{h}{\rho} \le \min_{\substack{0 \le s \le k \\ p+1 \le r \le m}} \left( \frac{\lambda_{s,p+1}}{\lambda_{s,r}} \right)^{1/(r-(p+1))} $$

where

$$ \lambda_{s,r} := \frac{r!}{(r-s)!} \sum_{l=\max\{s-r+p+1,0\}}^{s} \binom{s}{l} \binom{r-s}{p+1-l} $$

Then,

$$ \|u - (u^h + u^\rho)\|_{W_\infty^k} \le h^{p+1-k} [C_1 h^q + C_2 \rho^q] |u|_{W_\infty^{m+1}} \quad (19) $$

where $C_1$ and $C_2$ are independent of the finite element size, $h$, and the dilation parameter, $\rho$, of the meshless approximation.

The standard definition of the norm $\| \cdot \|_{W_\infty^k}$ is used, namely

$$ \|u\|_{W_\infty^k} = \sum_{s=0}^{k} |u|_{W_\infty^s} = \sum_{s=0}^{k} \sum_{|\alpha|=s} \max_{x \in \tilde{\Omega}} \left| \frac{\partial^{|\alpha|} u}{\partial x_1^{\alpha_1} \partial x_2^{\alpha_2} \dots \partial x_n^{\alpha_n}} \right| $$

# 6. NUMERICAL EXAMPLES

## 6.1. Coupled EFG-FEM

In this section a coupled FE-EFG approximation is employed with a simple example, the interpolation of $u(x) = \sin(\pi x)$ in $\Omega = [-1, 1]$. Linear elements are employed ($p = 1$) and the nodes in $\Omega^\rho = [-1, 0)$ are replaced by particles. Consistency of order one is enforced everywhere.

Figure 14 shows, on the left, the interpolation functions. The shape functions, $N^\rho$, associated to particles, denoted by asterisks, are plotted with a solid line. The finite element interpolation functions, $N^h$, are depicted with dashed lines and the position of the nodes by circles. The transition region $\tilde{\Omega}$ is $[-0.25, 0]$. Figure 14 also shows, on the right, the result of such an interpolation. The approximation $u^h + u^\rho$ is plotted with a solid line and the error, $u - (u^h + u^\rho)$, with a dashed line. It is important to notice the special profile adopted by the shape function associated to the first particle (particle at $x = -0.25$): on the left it is similar to the particle positioned at the boundary of the domain, while in $\tilde{\Omega}$ it looks like a standard linear finite element interpolation function.

In this case the approximation is similar to the one proposed in Reference [5]. However, here there is no need to define any ramp function. Moreover, the same formulation can be employed with a particles distribution such that the transition region $\tilde{\Omega}$ includes more than one particle. For instance, Figure 15 shows both the shape functions and the interpolation with its associated error for a different distribution of particles. In particular, now the transition region includes three particles (one on its boundary and two in the interior of $\tilde{\Omega}$). The larger number of particles (with their associated smaller dilation parameter) induces a better approximation in $\Omega^\rho = [-1, 0]$.
<figure><img src="image_16.png" /><figcaption>Figure 14. Approximation functions—4 particles and 5 nodes—(left) and interpolation result, $u^\rho + u^h$, with error distribution (right).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15. Approximation functions—12 particles and 5 nodes—(left) and interpolation result, $u^\rho + u^h$, with error distribution (right).</figcaption></figure>

## 6.2. Coupled and enriched EFG-FEM

Coupling and enrichment can be employed together. In this case, particles are added and element removed without any particular restriction. Function $u(x) = \sin(2\pi x)$ in $\bar{\Omega} = [-1, 1]$ is interpolated. As shown in Figure 16 four different regions are present: in $[-1, -0.5]$ only particles have an influence, in $[-0.5, 0]$ particles and a non-complete basis of finite elements are present, in $[0, 0.5]$ both particles and complete finite elements are used, finally, in $[0.5, 1]$ only finite elements have an influence. Consistency is not uniform in this case, in $\bar{\Omega} \setminus \Omega^\rho = [0.5, 1]$ the finite element interpolation controls the order of consistency, $m = p = 1$. But in the meshless area of influence, i.e. $\bar{\Omega}^\rho = [-1, 0.5]$, the order of consistency required is $m = 2$.

Figure 17 shows the interpolation results obtained with the particle distribution of Figure 16. Six particles and five nodes have been used, their associated shape functions are shown in
<figure><img src="image_18.png" /><figcaption>Figure 16. Approximation functions: 6 particles and 5 nodes.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17. Mixed interpolation with 6 particles and 5 nodes.</figcaption></figure>

<table><caption>Table I. Measures of error for 6 particles and 5 nodes.</caption><thead><tr><th></th><th>[-1, -0.5]</th><th>[-0.5, 0]</th><th>[0, 0.5]</th><th>[0.5, 1]</th></tr></thead><tbody><tr><td>Error in L<sub>2</sub> norm</td><td>0.059</td><td>0.098</td><td>0.073</td><td>0.107</td></tr><tr><td>Error in L<sub>∞</sub> norm</td><td>0.124</td><td>0.194</td><td>0.160</td><td>0.209</td></tr><tr><td>DOF EFG+FEM</td><td>2+0</td><td>2+1</td><td>2+3</td><td>0+3</td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>Figure 18. Mixed interpolation with 11 particles and 5 nodes.</figcaption></figure>

<table><caption>Table II. Measures of error for 11 particles and 5 nodes.</caption><thead><tr><th></th><th>[-1, -0.5]</th><th>[-0.5, 0]</th><th>[0, 0.5]</th><th>[0.5, 1]</th></tr></thead><tbody><tr><td>Error in L<sub>2</sub> norm</td><td>0.015</td><td>0.027</td><td>0.036</td><td>0.107</td></tr><tr><td>Error in L<sub>∞</sub> norm</td><td>0.048</td><td>0.052</td><td>0.088</td><td>0.209</td></tr><tr><td>DOF EFG+MEF</td><td>4+0</td><td>3+1</td><td>3+3</td><td>0+3</td></tr></tbody></table>

Figure 16. It is important to note that, as expected, the interpolation functions are hierarchical. The error in each region can be found in Table I.

In [0.5,1] linear finite elements induce the larger error. In [0,0.5] the error is reduced with an 'h-p' refinement: particles are added and the order of consistency is increased. In both regions the finite element interpolation is complete. The price, in the number of degrees of freedom, is considerable. Similar results are obtained if the number of particles is increased. Figure 18 shows the results with the same finite element mesh, the same orders of consistency and 11 particles. The dilation parameter $\rho$ is reduced by a half. The error measures can be found in Table II.

This example also shows the influence of a coarse finite element mesh when the number of particles is increased. This point was already discussed in the error analysis. In [-1, -0.5] and in [0, 0.5] the distribution of particles is similar. In the former the precision is higher albeit that the number of degrees of freedom is lower than in the other region. In the latter the complete finite element interpolation introduces extra degrees of freedom but the error does not decrease. As previously noted, see Section 5, if the finite element mesh is too coarse an increase in the number of particles does not reduce the error.

If the finite element mesh is enriched with meshless approximations, the coefficients associated to the finite element shape functions maintain their physical meaning. The meshless shape functions are hierarchical. However, convergence can only be achieved on a coarse mesh if the order of consistency is increased, i.e. adding more particles without any increase in $m$ does not suffice.
### 6.3. Finite element enrichment with EFG in a 2D scalar problem

The Poisson equation with Neumann and Dirichlet boundary conditions is solved next. The problem statement is

$$ 
\begin{aligned}
\Delta u &= -f, & \Omega &= (0,1) \times (0,1) \\
\nabla u \cdot \mathbf{n} &= q_0, & \Gamma_n &= \bar{\Omega} \cap \{y=0\} \\
u &= u_0, & \Gamma_d &= \partial\Omega \setminus \Gamma_n
\end{aligned}
$$

where **n** is the outward unit normal vector. The source term, $f$, and the boundary conditions, $q_0$ and $u_0$, are chosen such that $u(x) = e^{-(6(x+y-1))^2}$ is the solution. Plate 1 depicts this solution (left) and a cross-section on the plane $y = x$. Essential boundary conditions are imposed using Lagrange multipliers which are interpolated using the $C^0$ finite element interpolation functions along the boundary.

Plates 2 and 3 show the finite element mesh, the solution and the error distribution. An 8 × 8 quadrilateral mesh with bilinear finite elements (Q1) has been used. The error is larger along the diagonal $x + y = 1$ and the error measure in the maximum norm ($L_\infty$ norm) is 0.1707.

In order to improve the approximation, the finite element mesh is enriched adding particles and imposing an order of consistency $m = 2$. Plate 4 shows the finite element mesh and the distributions of particles. The error of the mixed approximation is also plotted in the same figure and with the same scale used in Plate 2. In fact the measure in the maximum norm is now: 0.0204.

Finally, Plate 5 presents the mixed approximation. The finite element approximation, $u^h$ (top), is improved by a particle contribution, $u^\rho$ (centre), which induces the final mixed approximation, $u^h + u^\rho$ (bottom).

### 6.4. Finite element enrichment with EFG in non-linear computational mechanics

This example reproduces the finite element enrichment with EFG in a nonlinear computational problem. A rectangular specimen with an imperfection is loaded, see References [19, 20]. It has two axes of symmetry, a bilinear elastoplastic material is considered, and plane strain conditions are assumed. Figure 19 presents the problem statement with the material properties.

This problem has been solved with standard eight-noded quadrilateral elements. Moreover, an adaptive error analysis [20, 21] has been conducted up to convergence. The final mesh and its equivalent inelastic strain distribution is shown in Plate 6 (left). This mesh has 2022 d.o.f. and a relative error (measured in energy norm) of 0.18 per cent.

The same example has also been solved with element-free Galerkin. In order to obtain comparable results, the distribution of particles coincides with the distribution of nodes in the previous finite element mesh; and consistency of order two is required. Thus, the number of degrees of freedom (d.o.f.) is also 2022. Plate 6 (right) shows the distribution of particles and inelastic strains.

Results degrade drastically if a coarse mesh of quadrilateral bilinear finite elements (308 d.o.f.) is employed, see Plate 7. However, when particles are added (308+906=1214 d.o.f.) and the order of consistency is increased ($m = 2$), the correct distribution of inelastic strains is recovered, see Plate 7. Note that, the final finite element mesh in Plate 6 (left) was obtained after an iterative process which needed for each iteration the generation of a new mesh. In this final example, Plate 7, the original mesh is maintained and particles are added where they are needed.

Finally, Figure 20 shows the evolution of the inelastic strains along the direction (A-A') for every configuration studied. Section (A-A') is plotted in Figure 19.
<figure><img src="image_21.png" /><figcaption>Figure 19. Problem statement: rectangular specimen with one centred imperfection.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 20. Force versus displacement (left) and evolution of the equivalent inelastic strain along (A-A') for each approximation (right).</figcaption></figure>

## 7. CONCLUSIONS

This paper develops a mixed interpolation, it is based on finite element and meshless methods. In fact, it is an extension of previous published papers by Liu, Belytschko and coworkers [5, 14, 16] with a unified formulation generalizable to any spatial order ($p$ or $m$) and with its corresponding convergence analysis. Two different cases have been studied: coupled finite elements with EFG, $p=m$, or enrichment of finite elements with EFG, $m>p$. For the sake of clarity, EFG has been used as the meshless method. However, generalization to RKPM is straightforward.
<figure><img src="image_23.png" /><figcaption>Plate 1. Analytical solution and section along $y = x$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Plate 2. Finite element mesh and error distribution.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Plate 3. Approximation with 8 x 8 Q1 finite elements.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Plate 4. Finite element mesh enriched with particles and error distribution of the mixed approximation.</figcaption></figure>
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /><figcaption>Plate 5. Finite element contribution $u^h$, enrichment with EFG $u^p$ and mixed approximation $u^h + u^p$.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Plate 6. Final mesh with its corresponding equivalent inelastic strain for a standard finite element (8 noded elements) computation (left) and distribution of particles with its inelastic strain distribution for EFG (right).</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Plate 7. Coarse finite element mesh (Q1 elements) with its corresponding equivalent inelastic strain (left) and mixed interpolation with its equivalent inelastic strain distribution (right).</figcaption></figure>
The first case allows to implement Dirichlet boundary conditions in a standard finite element context. In fact, this was proposed by Belytschko et al. [5] and here a simple generalization avoids the use of any ramp function and the need of substituting nodes per particles. That is, particles can be added arbitrarily in the region of the computational domain where the finite element interpolation is not complete. This ensures continuity of the solution (no coupling via Lagrange multipliers is imposed) and also enforces a uniform order of consistency (and thus of convergence) everywhere in the computational domain. The convergence properties of the mixed approximation are similar to those of the finite element method or element-free Galerkin.

The second case, enrichment of finite elements with EFG, allows to improve the accuracy of finite elements where needed in an adaptive process without any remeshing. The a priori error bounds of this formulation are illustrated with numerical examples and stated in a formal setting. They indicate that both the element size, *h*, and the dilation parameter, $\rho$, influence the convergence, as well as the order of the finite element interpolation, *p*, and the increase of order of consistency, *q*, due to the added particles. Moreover, the a priori bound shows that *h*, $h-p$ or $\rho-q$ refinements will induce convergence but $\rho$ refinement on its own will fail to decrease arbitrarily the error of the approximation. That is, convergence cannot be attained by simply adding particles and thus reducing the dilation parameter, an increase in the order of consistency is needed.

Both cases are illustrated with several examples. They show the applicability of the proposed formulation in standard linear and non-linear boundary value problems.

## REFERENCES

1. Liu WK, Chen Y. Wavelet and multiple scale reproducing kernel methods. *International Journal for Numerical Methods in Fluids* 1995; **21**:901-931.

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

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

4. Belytschko T, Organ D. Element-free Galerkin methods for dynamic fracture in concrete. In *Computer Plasticity. Fundamentals and Applications*, Owen DRJ, Oñate E, Hilton E (eds), 1997; 304-321.

5. Belytschko T, Organ D, Krongauz Y. A coupled finite element-free Galerkin method. *Computational Mechanics* 1995; **17**:186-195.

6. Belytschko T, Tabbara M. Dynamic fracture using element-free galerkin methods. *International Journal for Numerical Methods in Engineering* 1996; **39**:923-938.

7. Lu YY, Belytschko T, Gu L. A new implementation of the element free Galerkin method. *Computer Methods in Applied Mechanics and Engineering* 1994; **113**:397-414.

8. Organ D, Fleming M, Terry T, Belytschko T. Continuous meshless approximations for nonconvex bodies by diffraction and transparency. *Computational Mechanics* 1996; **18**:225-235.

9. Kulasegaram S, Bonet J. Corrected smooth particle hydrodynamics method for metal forming simulations. In *Simulation of Materials Processing: Theory, Methods and Applications*, Huétink J, Baaijens FPT (eds), 1998; 137-142.

10. Vila JP. On particle weighted methods and smooth particle hydrodynamics. *Mathematical Models and Methods in Applied Sciences* 1999; **9**:161-209.

11. Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. *Computer Methods in Applied Mechanics and Engineering* 1996; **139**:3-47.

12. Liu WK, Chen Y, Jun S, Chen JS, Belytschko T, Pan C, Uras RA, Chang CT. Overview and applications of the reproducing Kernel particle methods. *Archives of Computational Methods in Engineering, State of the Art Reviews* 1996; **3**:3-80.

13. Hegen D. Element free Galerkin methods in combination with finite element approaches. *Computer Methods in Applied Mechanics and Engineering* 1996; **135**:143-166.

14. Liu WK, Uras RA, Chen Y. Enrichment of the finite element method with reproducing Kernel particle method. *Journal of Applied Mechanics* 1997; **64**:861-870.

15. Liu WK, Belytschko T, Oden JT (eds), Meshless methods. *Computer Methods in Applied Mechanics and Engineering* 1996; **139**:1-440.

16. Liu WK, Li S, Belytschko T. Moving least square reproducing kernel methods. (I) Methodology and convergence. *Computer Methods in Applied Mechanics and Engineering* 1997; **143**:113-154.
17. Huerta A, Fernández-Méndez S, Díez P. Enrichissement des interpolations d'éléments finis en utilisant des méthodes de particules. *Annales de la Faculté des Sciences de Toulouse* (in French), to appear.

18. Fernández-Méndez S, Díez P, Huerta A. Convergence of finite elements enriched with meshless methods. *SIAM Journal of Numerical Analysis*, submitted.

19. Díez P, Arroyo M, Huerta A. Adaptivity based on error estimation for viscoplastic softening materials. *Mechanics of Cohesive-Frictional Materials* 2000; **5**:87–112.

20. Huerta A, Díez P. Error estimation including pollution assessment for nonlinear finite element analysis. *Computer Methods in Applied Mechanics and Engineering* 2000; **181**:21–41.

21. Huerta A, Rodríguez-Ferran A, Díez P, Sarrate J. Adaptive finite element strategies based on error analysis. *International Journal for Numerical Methods in Engineering* 1999; **46**:1803–1818.