# Enriched Reproducing Kernel Approximation: Reproducing Functions with Discontinuous Derivatives

Pierre Joyot, Jean Trunzler, Francisco Chinesta

► To cite this version:

Pierre Joyot, Jean Trunzler, Francisco Chinesta. Enriched Reproducing Kernel Approximation: Reproducing Functions with Discontinuous Derivatives. M. Griebel; M. Schweitzer. Meshfree Methods for Partial Differential Equations II. Lecture Notes in Computational Science and Engineering, 43, 2005, 10.1007/3-540-27099-X_6. hal-01633391

HAL Id: hal-01633391

https://hal.science/hal-01633391

Submitted on 12 Nov 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Enriched Reproducing Kernel Approximation: Reproducing Functions with Discontinuous Derivatives

Pierre Joyot<fnref n="1" />*, Jean Trunzler<fnref n="1" />,<fnref n="2" />**, and Fransisco Chinesta<fnref n="2" />***

**Abstract** In this paper we propose a new approximation technique within the context of meshless methods able to reproduce functions with discontinuous derivatives. This approach involves some concepts of the reproducing kernel particle method (RKPM), which have been extended in order to reproduce functions with discontinuous derivatives. This strategy will be referred as Enriched Reproducing Kernel Particle Approximation (E-RKPA). The accuracy of the proposed technique will be compared with standard RKP approximations (which only reproduces polynomials).

## 1 Introduction

Meshless methods are an appealing choice to develop functional approximations (with different degrees of consistency and continuity) without a mesh support. Thus, this kind of techniques seem to be specially appropriated for treating 3D problems involving large displacements, due to the fact that the approximation is constructed only from the cloud of nodes whose position evolve during the material deformation. In this manner neither remeshing nor fields projections are a priori required.

In order to adapt the approximation for introducing some known information associated to the searched solution different possibilities exist. The first one lies in the enrichment of standard approximations in the framework of the partition of unity (PU), originally proposed by Babuška and Melenk [3], and which is at the origin of the partition of unity finite elements (PUFEM) [12]. Generalized finite elements (G-FEM) [15] or extended finite elements (X-FEM) [13, 16, 17, 9] are two members of this family. In the extended finite elements the approximation is locally enhanced by introducing appropriate functions describing the known behavior of the problem solution.

<footnote n="1">¹ LIPSI-ESTIA, Technopole Izarbel, 64210 Bidart, France.</footnote>

<footnote n="2">² LMSP, 151 Bd. de l'Hôpital, 75013 Paris, France.</footnote>

<footnote>* p.joyot@estia.fr</footnote>

<footnote>** j.trunzler@estia.fr</footnote>

<footnote>*** francisco.chinesta@paris.ensam.fr</footnote>
The second possibility lies in the introduction of some information related to the problem solution within the approximation functional basis. For this purpose, different reproduction conditions are enforced in the construction of the approximation functions. This approach has been widely used in the context of the moving least squares approximations used in the diffuse meshless techniques [2] as well as in the element free Galerkin method [4]. Thus, all the known information related to the problem solution can be introduced in the functional approximation. Very accurate results were obtained for example in fracture mechanics by introducing the crack tip behavior into the approximation basis [6].

In this work we propose a numerical strategy, based on the reproducing kernel particle techniques, able to construct approximation functions with discontinuous derivatives on fixed or moving interfaces. This problem was treated in the context of the partition of unity by Kronggauz et al. [10]. In our approach the size of the discrete system of equations remains unchanged because no additional degrees of freedom are introduced related to the enrichment. However, the fact of enriching the approximation implies, as shown later, a bigger moment matrix, with worse properties, but in some cases this enrichment is only local, and in any case, the moment matrix has a low dimension.

The starting point of our development is the reproducing kernel particle approximation (RKPA). The RKP approximation was introduced by Liu et al. [11] for enforcing some degree of consistency to standard smooth particle approximations, i.e. they proved that starting from a SPH (smooth particle hydrodynamics) approximation [7] it is possible to enhance the kernel function for reproducing a certain degree of polynomials. We have extended or generalized this procedure in order to reproduce any function, and more concretely, functions involving discontinuous derivatives.

Moreover, the strategy here described, can be easily coupled with a level set description of the interface location [14] (recently introduced in the framework of the extended finite elements [17]), which allows to capture its position when it evolves in time. Thus, both the enrichment and the interface tracking are made in a transparent way for the user.

This paper consists of 6 sections, in which we summarize the reproducing kernel approximation (section 2) which allows to treat in section 3 the reproduction of a function with discontinuous derivatives. In sections 4 and 5 we describe and illustrate the properties of the resulting approximation. Finally, in section 6 we focus on a test problem which allows to conclude about the accuracy and potentiality of the proposed strategy.
## 2 Enriched Reproducing Kernel Particle Approximation

Let $\Omega$ be the domain where the problem is defined. The points within this domain will be noted by $x$ or $s$. For the sake of simplicity, from now on, we will focus in the 1D case, but all the results have a direct 2D or 3D counterpart.

### 2.1 Reproduction Conditions

The approximation $f^a(x)$ of $f(x)$ is built from the integral convolution integral

$$f^a(x) = \int_{\Omega} \Phi(x-s,h)f(s)d\Omega, \quad (2.1)$$

where $\Phi(x-s,h)$ is the kernel function and $h$ a parameter defining the size of the approximation support.

The main aim of this work is to enforce the reproduction of a general function that we can write in the form of a polynomial plus another function noted by $f^e(x)$:

$$f^a(x) = a_0 + a_1 x + \dots + a_n x^n + f^e(x). \quad (2.2)$$

In the following paragraphs we analyze the required properties of the kernel function $\Phi(x-s,h)$ for reproducing a function expressed by (2.2). From Eq. (2.1), the reproduction of a constant function $a_0$ is given by

$$\int_{\Omega} \Phi(x-s,h)a_0 d\Omega = a_0, \quad (2.3)$$

which implies

$$\int_{\Omega} \Phi(x-s,h)d\Omega = 1, \quad (2.4)$$

which constitutes the partition of unity property.

Now, the required condition to reproduce a linear function $f^a(x) = a_0 + a_1 x$ is

$$\int_{\Omega} \Phi(x-s,h)(a_0 + a_1 s) d\Omega = a_0 + a_1 x. \quad (2.5)$$

By using the partition of unity (2.4), Eq. (2.5) can be rewritten as

$$\left\{ \begin{array}{l} \int_{\Omega} \Phi(x-s,h) d\Omega = 1 \\ \int_{\Omega} \Phi(x-s,h) s d\Omega = x, \end{array} \right. \quad (2.6)$$

which implies the linear consistency of the approximation. Repeating this reasoning, we can write the $n$-order consistency as

$$\left\{ \begin{array}{l} \int_{\Omega} \Phi(x-s,h) d\Omega = 1 \\ \int_{\Omega} \Phi(x-s,h) s d\Omega = x \\ \vdots \\ \int_{\Omega} \Phi(x-s,h) s^n d\Omega = x^n \end{array} \right. \quad (2.7)$$

and consequently, the reproduction of the function given by (2.2) implies

$$\int_{\Omega} \Phi(x-s,h)(a_0 + a_1 s + \dots + a_n s^n + f^e(s))d\Omega =$$
$$= a_0 + a_1 x + \dots + a_n x^n + f^e(x), \qquad (2.8)$$

from which it results

$$\left\{ \begin{array}{l} \displaystyle \int_{\Omega} \Phi(x-s,h) d\Omega = 1 \\ \displaystyle \int_{\Omega} \Phi(x-s,h) s d\Omega = x \\ \vdots \\ \displaystyle \int_{\Omega} \Phi(x-s,h) s^n d\Omega = x^n \\ \displaystyle \int_{\Omega} \Phi(x-s,h) f^e(s) d\Omega = f^e(x) \end{array} \right. . \qquad (2.9)$$

Chen et al. [5] applied the procedure proposed by Liu et al. [11] for enforcing n-order consistency, however their procedure does not allow to enforce directly the reproduction condition associated with $f^e(x)$.

## 2.2 The Moment Matrix

We will note by $f^r(x)$ the approximation function verifying the conditions (2.9). Usually a cubic spline is considered as kernel function, and consequently the conditions given by Eq. (2.9) are not satisfied. Liu et al. [11] propose the introduction of a correction function $C(x, x-s)$ for satisfying the reproduction conditions. In our case we consider the more general form $C(x, s, x-s)$ whose pertinence will be discussed later. Thus $f^r(x)$ will be expressed by

$$f^r(x) = \int_{\Omega} C(x, s, x-s) \Phi(x-s, h) f(s) d\Omega, \qquad (2.10)$$

where $C(x, s, x-s)$ is assumed to have the following form

$$C(x, s, x-s) = \mathbf{H}^T(x, s, x-s)\mathbf{b}(x), \qquad (2.11)$$

where $\mathbf{H}^T(x, s, x-s)$ represents the vector containing the functions considered in the approximation basis, and $\mathbf{b}(x)$ is a vector containing unknown functions that will be determined for satisfying the reproduction conditions. Thus, Eq. (2.9) can be rewritten as

$$\left\{ \begin{array}{l} \displaystyle \int_{\Omega} \mathbf{H}^T(x, s, x-s) \mathbf{b}(x) \Phi(x-s, h) d\Omega = 1 \\ \displaystyle \int_{\Omega} \mathbf{H}^T(x, s, x-s) \mathbf{b}(x) \Phi(x-s, h) s d\Omega = x \\ \vdots \\ \displaystyle \int_{\Omega} \mathbf{H}^T(x, s, x-s) \mathbf{b}(x) \Phi(x-s, h) s^n d\Omega = x^n \\ \displaystyle \int_{\Omega} \mathbf{H}^T(x, s, x-s) \mathbf{b}(x) \Phi(x-s, h) f^e(s) d\Omega = f^e(x) \end{array} \right. . \qquad (2.12)$$

*Remark:* Usual RKP approximations consider $\mathbf{H}^T(x-s)$, but we retain the more general form $\mathbf{H}^T(x, s, x-s)$ in order to reproduce more general functions.

In fact, the reproduction conditions must be enforced in a discrete form. For this purpose we consider NP points (also refereed as nodes) which allow to compute the discrete form of Eq. (2.12), i.e.
$$ egin{cases} 
\sum_{I=1}^{NP} \mathbf{H}^T(x, x_I, x-x_I) \mathbf{b}(x) \Phi(x-x_I, h) \Delta x_I = 1 \\
\sum_{I=1}^{NP} \mathbf{H}^T(x, x_I, x-x_I) \mathbf{b}(x) \Phi(x-x_I, h) x_I \Delta x_I = x \\
\vdots \\
\sum_{I=1}^{NP} \mathbf{H}^T(x, x_I, x-x_I) \mathbf{b}(x) \Phi(x-x_I, h) x_I^\eta \Delta x_I = x^\eta \\
\sum_{I=1}^{NP} \mathbf{H}^T(x, x_I, x-x_I) \mathbf{b}(x) \Phi(x-x_I, h) f^e(x_I) \Delta x_I = f^e(x)
\end{cases} 
, \quad (2.13) $$

that in a matrix form results

$$ \left[ \sum_{I=1}^{NP} \mathbf{R}(x_I) \mathbf{H}^T(x, x_I, x-x_I) \Phi(x-x_I, h) \Delta x_I \right] \mathbf{b}(x) = \mathbf{R}(x), \quad (2.14) $$

where $\mathbf{R}(x)$ is the reproduction vector

$$ \mathbf{R}^T(x) = [1, x, \dots, x^n, f^e(x)]. \quad (2.15) $$

Eq. (2.14) allows the computation of vector $\mathbf{b}(x)$,

$$ \mathbf{b}(x) = \mathbf{M}(x)^{-1} \mathbf{R}(x), \quad (2.16) $$

where the moment matrix $\mathbf{M}(x)$ is defined by

$$ \mathbf{M}(x) = \sum_{I=1}^{NP} \mathbf{R}(x_I) \mathbf{H}^T(x, x_I, x-x_I) \Phi(x-x_I, h) \Delta x_I. \quad (2.17) $$

This moment matrix differs of the usual moment matrix proposed in [11], and in fact it becomes non symmetric.

## 2.3 Discrete Form of the Approximation Function

The discrete form $f^r(x)$ of $f^a(x)$ derives from Eqs. (2.10), (2.11) and (2.16)

$$ 
\begin{aligned}
f^r(x) &\approx \sum_{I=1}^{NP} \mathbf{H}^T(x, x_I, x-x_I) \mathbf{M}(x)^{-1} \mathbf{R}(x) \Phi(x-x_I, h) f(x_I) \Delta x_I \\
&= \sum_{I=1}^{NP} \psi_I(x) f_I,
\end{aligned} 
\quad (2.18) $$

where $\psi_I$ is the enriched RKP approximation shape function

$$ \psi_I(x) = \mathbf{H}^T(x, x_I, x-x_I) \mathbf{M}(x)^{-1} \mathbf{R}(x) \Phi(x-x_I, h) \Delta x_I. \quad (2.19) $$

As in the classical RKPM we take $\Delta x_I = 1$. Different quadrature rules exist and they have been tested in [1] without a significant incidence on the reproducing condition accuracy.

## 3 The Case of a Function with Discontinuous Derivatives

Let $\Omega$ be the domain where the problem is defined and $\Gamma_d$ a point, curve or surface (in 1D, 2D and 3D respectively) where the normal function derivative becomes discontinuous. We assume that this discontinuity curve splits the domain in two subdomains $\Omega_0$ and $\Omega_1$ (see Figure 3.1)
<figure><img src="image_2.png" /><figcaption>**Figure 3.1.** Problem domain containing an interface with a discontinuous normal derivative.</figcaption></figure>

$$ \begin{aligned} \Omega_0 \cup \Omega_1 \cup \Gamma_d &= \Omega, \\ \Omega_0 \cap \Omega_1 &= \emptyset, \end{aligned} $$

where $\Gamma_d$ represents the curve along which the normal derivative of some approximation functions become discontinuous.

The enrichment function $f^e(\boldsymbol{x})$ which will be introduced in the reproduction vector $\boldsymbol{R}(\boldsymbol{x})$ must satisfy the following conditions

$$ \begin{aligned} f^e(\boldsymbol{x}) &\in \mathcal{C}^0(\boldsymbol{x} \in \Omega), \\ f^e(\boldsymbol{x}) &\in \mathcal{C}^1(\boldsymbol{x} \in \Omega_0), \\ f^e(\boldsymbol{x}) &\in \mathcal{C}^1(\boldsymbol{x} \in \Omega_1), \\ \mathrm{Grad} \, f^e(\boldsymbol{x}) \cdot \boldsymbol{n} &\in \mathcal{C}^{-1}(\boldsymbol{x} \in \Gamma_d), \end{aligned} \quad (3.1) $$

where $\boldsymbol{n}$ denotes the unit outward vector defined on the curve $\Gamma_d$.

To locate the discontinuity curve $\Gamma_d$ we have use of a level set function $\Theta(\boldsymbol{x})$ defined as the signed distance from $\boldsymbol{x}$ to the interface $\Gamma_d$. Thus

$$ \Theta(\boldsymbol{x}) = \begin{cases} \Theta(\boldsymbol{x}) < 0 & \text{if } \boldsymbol{x} \in \Omega_0 \\ \Theta(\boldsymbol{x}) > 0 & \text{if } \boldsymbol{x} \in \Omega_1 \\ \Theta(\boldsymbol{x}) = 0 & \text{if } \boldsymbol{x} \in \Gamma_d \end{cases}. \quad (3.2) $$

Thus, for satisfying the conditions (3.1) the enrichment function $f^e(\boldsymbol{x})$ could be assumed in the form

$$ f^e(\boldsymbol{x}) = H_0(\Theta(\boldsymbol{x}))\Theta(\boldsymbol{x}), \quad (3.3) $$

where $H_0(\boldsymbol{x})$ represents the usual Heaviside function

$$ \left\{ \begin{array}{ll} H_0(\Theta(\boldsymbol{x})) = 1 & \text{if } \Theta(\boldsymbol{x}) \ge 0 \\ H_0(\Theta(\boldsymbol{x})) = 0 & \text{if } \Theta(\boldsymbol{x}) < 0 \end{array} \right. . \quad (3.4) $$

# 4 Properties of the Moment Matrix

If the moment matrix is computed as described in section 2 everywhere, it can become singular, as we will prove later. We can define two boundaries in the discontinuity curve neighborhood $\Gamma_0$ and $\Gamma_1$ (as illustrated in figure 4.2)
<figure><img src="image_3.png" /><figcaption>Figure 4.2. Enriched and non-enriched domains.</figcaption></figure>

such that now the domain $\Omega$ consists of four subdomains $\Omega_{c0}$, $\Omega_{e0}$, $\Omega_{c1}$ and $\Omega_{e1}$, where classical and enriched reproducing kernel approximation will be defined with continuity along $\Gamma_0$ and $\Gamma_1$ as proved later. Thus, we start from the relations

$$ \begin{aligned} \Omega_0 &= \Omega_{c0} \cup \Omega_{e0} \cup \Gamma_0, \\ \Omega_1 &= \Omega_{c1} \cup \Omega_{e1} \cup \Gamma_1. \end{aligned} $$

Obviously, $\Gamma_0$ and $\Gamma_1$ define the points, curves or surfaces (1D, 2D and 3D respectively) in $\Omega_0$ and $\Omega_1$ respectively, through which the determinant of the moment matrix changes from zero to a non-null value.

We define $\mathbf{R}_c^T(x) = [1, x, \dots, x^n]$, $\mathbf{R}^T(x) = [\mathbf{R}_c^T(x), H_0(\Theta(x))\Theta(x)]$, and in a similar manner $\mathbf{H}_c^T(x-x_I) = [1, x-x_I, \dots, (x-x_I)^n]$, and $\mathbf{H}^T(x, x_I, x-x_I) = [\mathbf{H}_c^T(x-x_I), H_0(\Theta(x_I))\Theta(x_I)]$.

The following properties of the moment matrix can be deduced:

*Property 1:* The inverse of the moment matrix exists in the domains $\Omega_{e0}$ and $\Omega_{e1}$, i.e. in the subdomains whose boundary involves the discontinuity curve $\Gamma_d$. Then, the moment matrix results

$$ M(x) = \sum_{I=1}^{NP} \begin{bmatrix} \mathbf{R}_c(x_I)\mathbf{H}_c^T(x-x_I) & \mathbf{R}_c(x_I)H_0(\Theta(x_I))\Theta(x_I) \\ H_0(\Theta(x_I))\Theta(x_I)\mathbf{H}_c^T(x-x_I) & H_0(\Theta(x_I))^2\Theta(x_I)^2 \end{bmatrix} \Phi(x-x_I,h)\Delta x_I. $$

*Property 2:* In the domain $\Omega_{c0}$, $H_0(\Theta(x_I)) = 0$, $\forall I$. Thus, the moment matrix becomes singular

$$ M(x) = \sum_{I=1}^{NP} \begin{bmatrix} \mathbf{R}_c(x_I)\mathbf{H}_c^T(x-x_I) & 0 \\ 0 & 0 \end{bmatrix} \Phi(x-x_I,h)\Delta x_I. $$

In this case we propose to replace this moment matrix for the one derived from a standard (non-enriched) approximation
$$ M_c(x) = \sum_{I=1}^{ND} R_c(x) H_c^T (x - x_I) \Phi(x - x_I, h) \Delta x_I. \quad (4.1) $$

*Property 3:* Within the domain $\Omega_{c1}$, $H_0(\Theta(x_I)) = 1$, $\forall I$. Thus, $M(x)$ results

$$ M(x) = \sum_{I=1}^{NP} \begin{bmatrix} R_c(x_I)H_c^T(x-x_I) & R_c(x_I)\Theta(x_I) \\ \Theta(x_I)H_c^T(x-x_I) & \Theta(x_I)^2 \end{bmatrix} \Phi(x-x_I, h) \Delta x_I. $$

The moment matrix is singular in 1D (as will be illustrated in section 5.1) as well as in some 2D and 3D cases. In this case we replace again the enriched moment matrix for the non-enriched one. Effectively, if the enriched moment matrix is not singular everywhere in $\Omega_{c1}$, then $\Gamma_1$ does not exist, i.e. $\Omega_{c1} = \emptyset$ et $\Omega_{e1} = \Omega_1$.

*Property 4:* When we approach the curve $\Gamma_0$ the shape functions computed from both the enriched and the non-enriched approximations are connected with continuity, i.e.

$$ \lim_{x \in \Omega_{c0} \to \Gamma_0} \psi_I(x) = \lim_{x \in \Omega_{e0} \to \Gamma_0} \psi_I(x), \quad \forall I. $$

And in a similar way when we approach the curve $\Gamma_1$ (when it exists)

$$ \lim_{x \in \Omega_{c1} \to \Gamma_1} \psi_I(x) = \lim_{x \in \Omega_{e1} \to \Gamma_1} \psi_I(x), \quad \forall I. $$

This property has important consequences because due to the continuity of the shape functions within the whole domain, the numerical integration of variational formulations does not require a specific treatment.

*Property 5:* Let $\zeta_I$ be the boundary of the support related to the node $x_I$. The following results can be stated:

The curves $\Gamma_0$ and $\Gamma_1$ belong to the union of the support boundary $\zeta_I$ of some nodes. In the 1D case $\Gamma_0$ and $\Gamma_1$ reduce to two points (as illustrated in section 5). In higher dimensions, both boundaries are composed by a series of arcs of circle (sphere in 3D) (when the support of the kernel function is assumed circular -spherical in 3D-).

*Property 6:* The boundary $\zeta_I$ associated with the node $x_I \in \Omega_1$ (respectively $x_I \in \Omega_0$) closest to the discontinuity interface $\Gamma_d$ defines the boundary $\Gamma_0$ (respectively $\Gamma_1$).

# 5 About the Resulting Shape Function and its Derivatives

In this section we focus in a one-dimensional problem in order to illustrate some of the previous properties. The signed distance is used for defining the level set function
$$ \Theta(x) = x - x_d, \qquad (5.1) $$

where $x_d$ represents the coordinate of the point where the discontinuity in the function derivative takes place.

## 5.1 Moment Matrix

For the sake of clarity, from now on, we restrict the consistency requirement to a first order. Thus, the moment matrix results

$$ M(x) = \sum_{I=1}^{NP} \tilde{M}_{I}(x) \Phi(x - x_{I}, h) $$

with

$$ \tilde{M}_{I}(x) = \begin{bmatrix} 1 & x-x_{I} & H_{0}(x_{I}-x_{d})(x_{I}-x_{d}) \\ x_{I} & x_{I}(x-x_{I}) & x_{I}H_{0}(x_{I}-x_{d})(x_{I}-x_{d}) \\ H_{0}(x_{I}-x_{d})(x_{I}-x_{d}) & H_{0}(x_{I}-x_{d})(x_{I}-x_{d})(x-x_{I}) & (H_{0}(x_{I}-x_{d})(x_{I}-x_{d}))^{2} \end{bmatrix}. $$

For illustrating the third property, we consider $x \in \Omega_{c1}$, which leads to the following expression of the moment matrix

$$ M(x) = \sum_{I=1}^{NP} \begin{bmatrix} 1 & x-x_{I} & (x_{I}-x_{d}) \\ x_{I} & x_{I}(x-x_{I}) & x_{I}(x_{I}-x_{d}) \\ (x_{I}-x_{d}) & (x_{I}-x_{d})(x-x_{I}) & (x_{I}-x_{d})^{2} \end{bmatrix} \Phi(x-x_{I},h). $$

We can notice that the third row is a linear combination of the two first. Thus, the determinant vanishes and the moment matrix becomes singular. This result does not depend on the degree of consistency assumed.

Now we are going to illustrate the position of the boundaries $\Gamma_0$ and $\Gamma_1$ through which the approximation changes from enriched to non-enriched.

In Figure 5.3 the vertical line indicates the position where the function derivative is expected to be discontinuous, and the horizontal one represents the enrichment domain $\Omega_{e1} \cup \Omega_{e0}$. The distance between two consecutive nodes (assumed regularly spaced) is $h = 1$ and the support radius is taken in the form $R = a \times h$ (where the value of $a$ is indicated in each subfigure).

It is well known that for a standard RKP approximation (without any enrichment), the linear consistency implies that each point must be contained by more than one nodal support (in other case the moment matrix becomes singular) [8]. In our case we need to add to this constraint the properties listed in the previous section.

From Figure 5.3(b) we can notice:

* Point A is in the domain $\Omega_{c0}$, i.e. $A \in \Omega_{c0}$ because any node whose support contains the point A is not in $\Omega_1$ (Property 2).

* $D \in \Omega_{c1}$ because any node whose support contains the point D is not in $\Omega_0$ (Property 3).
<figure><img src="image_4.png" /></figure>

(a) Node not located at the discontinuity: discontinuous $\Omega_{e0} \cup \Omega_{e1}$.

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

(b) Node not located at the discontinuity: continuous $\Omega_{e0} \cup \Omega_{e1}$.

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

(c) Node located at the discontinuity.

**Figure 5.3.** Enriched and non-enriched domains in a one-dimensional test problem.

• $B \in \Omega_{e0}$ and $C \in \Omega_{e1}$ (Property 1).

• The position of $\Gamma_1$ (respectively $\Gamma_0$) is defined by the support of the node within $\Omega_0$ (respectively $\Omega_1$) closest to the discontinuity interface $\Gamma_d$ (Properties 5-6).

If $a < 1.5$ there are points, candidates to belong to $\Omega_{e0} \cup \Omega_{e1}$, such that they belong only to the support of two nodes. In this way the approximation is only defined if $a > 1.5$ (Figure 5.3(a)).

When there is a node located at the discontinuity point this problem is avoided as depicted in Figure 5.3(c). Thus, the approximation is defined only if $a > 1$, like in the standard (non-enriched) RKP approximation.

## 5.2 Related Shape Functions

Figure 5.4 depicts the shape functions related to each node considered in the interval $\Omega = [-1, 1]$ with $x_d = 0$ and $NP = 10$ ($h = 2/9$). The support radius $a \times h$ is defined by the value of the parameter $a$, that has been fixed to $a = 2$. Thus, the boundaries defining the transition between the enriched and non-enriched domain are located at $x_0 = -1.5 \times h = -\frac{1}{3}$ and $x_1 = 1.5 \times h = \frac{1}{3}$. Figure 5.5 shows the associated derivatives, whereas Figure 5.6 depicts
<figure><img src="image_7.png" /><figcaption>**Figure 5.4.** Shape functions for a discontinuous derivative located at $x_d = 0$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.5.** Shape function derivatives for a discontinuous derivative located at $x_d=0$.</figcaption></figure>

similar results when $NP = 9$ and consequently a node is placed just at the discontinuity point.

*Remarks:*

* We have noticed that the shape functions are perfectly continuous through the interfaces between the enriched and non-enriched domains (*Property 4*). Moreover the shape function derivatives become discontinuous only along the interface $x_d = 0$.

* The shape functions are quite different when a node is or not located at the discontinuity point. In the first case (node located at $x_d$) the shape functions are very close to those obtained by Krongauz et al. [10] using a PU enrichment.
<figure><img src="image_9.png" /><figcaption>**Figure 5.6.** Shape functions for a discontinuous derivative located at $x_d = 0$ when a node is placed just at the discontinuity point.</figcaption></figure>

## 6 Numerical Example

In this section we consider a simple numerical test: the steady heat conduction problem (with a constant source term $g = 2$) defined in the interval $]-1, 1[$ which involves two homogeneous materials with different thermal conductivities

$$ \frac{d}{dx} \left( k(x) \frac{dT}{dx} \right) = g, \qquad (6.1) $$

where $T$ is the temperature field and $k$ the thermal conductivity defined as follows

$$ \begin{aligned} k(x) &= k_0 = 1 \quad \text{if } x \in \Omega_0 = ]-1, 0[, \\ k(x) &= k_1 = 10 \quad \text{if } x \in \Omega_1 = ]0, 1[. \end{aligned} \qquad (6.2) $$

The boundary conditions are defined by

$$ \left\{ \begin{array}{l} T(x = -1) = 0 \\ T(x = 1) = 1 \end{array} \right. . \qquad (6.3) $$

The problem is solved using the weak-form of Eq. (6.1), where the integrals were computed using a background mesh. Four integration cells are placed between two consecutive nodes, with 5 gauss integration points in each cell.

Computed temperatures with and without approximation enrichment as well as the exact temperature field are depicted in figure 6.7. In both cases (enriched and non-enriched) a first order consistency has been enforced, been the support radius $R = 2 \times h$. The related derivatives are shown in Figure 6.8. From both figures we can conclude about the better accuracy of the solution computed by using an enriched approximation, mainly in the temperature derivative in the discontinuity neighbourhood.
<figure><img src="image_10.png" /><figcaption>**Figure 6.7.** Exact and numerical temperature fields computed with and without approximation enrichment.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.8.** Exact and computed temperature derivatives.</figcaption></figure>

## 6.1 Convergence Analysis

In this section we compare the convergence of both enriched and non-enriched approximations, in the test problem just described. Two norms are considered, the $L_2$ and the $H_1$ norms. Moreover, the analysis has been performed for three different conductivity ratios $k_{01} = \frac{k_1}{k_0} = 10, 100, 1000$. Figure 6.9 illustrates these results.

In the case of non-enriched approximations (standard RKPA) we can notice from figure 6.9 that the order of the method is 1 and 0.5 (using the $L_2$ and $H_1$ norms respectively). Moreover, the error increases as the conductivity ratio increases.
On the contrary, when an enriched approximation is considered the method
results of order 2.4 and 1.5 (using the $L_2$ and $H_1$ norms respectively). More-
over, the error is independent of the conductivity ratio.

In conclusion, both the order of convergence and the solution error ob-
tained by using an enriched approximation are significantly better that the
ones obtained by using the standard non-enriched formulation.

<figure><img src="image_12.png" /><figcaption>Figure 6.9. Convergence analysis using the $L_2$ and the $H_1$ norms.</figcaption></figure>

7 Conclusion

In this work we have proposed a new strategy to enrich reproducing kernel
particle approximations in order to reproduce functions involving discontin-
uous derivatives through some interfaces. The computed results seems to be
very accurate (from the point of view of the order of convergence as well as
from the associated level of error) and it is very easy to implement.

In some cases the moment matrix becomes singular far from the discontinu-
ity location (in the order of the support size), and consequently the enrichment
becomes local in space. In any case, the size of the global discrete system is
not affected by the enrichment process, in contrast with the vast majority of
enrichment methods making use of the partition of unity.

The case of moving interfaces involving field discontinuities can be easily
treated coupling the technique described in this work with a level set descrip-
tion of the interface movement.

The extension to higher dimension problems involving complex disconti-
nuity geometries is a work in progress.
References

1. N. Aluru. A point collocation method based on reproducing kernel approxima-
tions. *Int. J. Numer. Meths. Eng.*, 47:1083–1121, 2000.

2. G. T. B. Nayrolles and P. Villon. Generalizing the finite element method: Dif-
fuse approximation and diffuse elements computational mechanics. *Journal of*
*Computational Mechanics*, 10(5):307–318, 1992.

3. I. Babuška and J. Melenk. The partition of unity method. *Int. J. Numer. Meths. Eng.*, 40:727–758, 1997.

4. T. Belytschko, Y. Kronggauz, D. Organ, and M. Fleming. Meshless methods: an overview and recent developments. *Comput. Meths. Appl. Mech. Engrg.*, 139:3–47, 1996.

5. J. Chen, C. Pan, C. Roque, and H.P.Wang. A lagrangian reproducing kernel particle method for metal forming. *Comput. Mech.*, 22:289–307, 1998.

6. M. Fleming, Y. Chu, B. Moran, and T. Belytschko. Enriched element-free galerkin methods for crack tip fields. *Int. J. Numer. Meths. Eng.*, 40, 1997.

7. R. Gingold and J. Monaghan. Smoothed particle hydrodynamics: Theory and application to non-spherical stars. *Monthly Notices Royal Astronomical Society*, 181:375–389, 1977.

8. W. Han and X. Meng. Error analysis of the reproducing kernel particle method. *Comput. Meths. Appl. Mech. Engrg.*, 190:6157–6181, 2001.

9. H. Ji, D. Chopp, and J. Dolbow. A hybrid extended finite/level set method for modeling phase transformations. *Int. J. Numer. Meths. Eng.*, 54:1209–1233, 2002.

10. Y. Kronggauz and T. Belytschko. Efg approximation with discontinuous derivatives. *Int. J. Numer. Meths. Eng.*, 41:1215–1233, 1998.

11. W. Liu, S. Jun, S. Li, J. Adee, and T. Belytschko. Reproducing kernel particle methods for structural dynamics. *Int. J. Numer. Meths. Eng.*, 38:1655–1679, 1995.

12. J. M. Melenk and I. Babuška. The partition of unity finite element method: Basic theory and applications. *Comput. Meths. Appl. Mech. Engrg.*, 139:289–314, 1996.

13. N. Moës and T. Belytschko. Extended finite element method for cohesive crack growth. *Engineering Fracture Mechanics*, 69:813–833, 2002.

14. J. Sethian. Evolution, implementation and application of level set and fast marching methods for advancing fronts. *J. Comput. Phys.*, 169:503–555, 2001.

15. T. Strouboulis, K. Corps, and I. Babuška. The generalized finite element method.
*Comput. Meths. Appl. Mech. Engrg.*, 190:4081–4193, 2001.

16. N. Sukumar, D. Chopp, N. Moës, and T. Belytschko. Modeling holes and in-
clusions by level sets in the extended finite-element method. *Comput. Meths.
Appl. Mech. Engrg.*, 190:6183–6200, 2001.

17. G. Ventura, J. Xu, and T. Belytschko. A vector level set method and new discontinuity approximations for crack growth by efg. *Int. J. Numer. Meths. Eng.*, 54:923–944, 2002.