# An Analysis of the Performance of Several Natural Neighbour-Based Galerkin Methods

Icíar Alfaro, Julien Yvonnet, Francisco Chinesta, Elías Cueto, Manuel Doblaré

► To cite this version:

Icíar Alfaro, Julien Yvonnet, Francisco Chinesta, Elías Cueto, Manuel Doblaré. An Analysis of the Performance of Several Natural Neighbour-Based Galerkin Methods. 9th ESAFORM Conference on Material Forming, 2006, Glasgow, United Kingdom. hal-00080987

HAL Id: hal-00080987

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

Submitted on 13 Mar 2018

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 CC BY 4.0 - Attribution - International License
# An analysis of the performance of several Natural Neighbour-based Galerkin methods

I. Alfaro¹, J. Yvonnet², F. Chinesta³, E. Cueto¹, M. Doblaré¹

¹Aragón Institute of Engineering Research. University of Zaragoza - María de Luna, s.n. 50018 Zaragoza, Spain

URL: http://mmcyte7.cps.unizar.es

e-mail: iciar@unizar.es; ecueto@unizar.es

²Laboratoire de Mécanique. Université Marne-la-Valée -5 Bd Descartes, 77454 Marne-la-vallée Cedex 2

URL: www.univ-mlv.fr

e-mail: jyvonnet@univ-mlv.fr;

³LMSP, ENSAM Paris - 151 Bvd. De l'Hôpital, 75013 Paris. France

URL: http://lms-web.paris.ensam.fr/lms

e-mail: francisco.chinesta@paris.ensam.fr

ABSTRACT: The Natural Element Method (NEM) has better accuracy than linear Finite Element Method when very distorted meshes are used. However the time needed to make the same calculation is higher in the meshless method than in FEM. Nevertheless, this computational cost includes, in the case of the meshless method, the user time employed in the mesh generation in the FEM. This article compares the cost in time using FEM and NEM with different shape functions and the accuracy of the obtained results.

Key words: NEM, meshless, Sibson, Laplace, Pseudo NEM, computing time.

## 1 INTRODUCTION

The natural element method (NEM) is a particular meshfree method which possesses advantages compared to the vast majority of other meshless approaches such as (a) Direct imposition of essential boundary conditions (b) Automatic adaptation of the support (domain of influence of a node) to the density and/or anisotropy of the surroundings cloud of nodes. Furthermore, it has been shown that the accuracy of the approximation is not much affected by the relative positions of the nodes. In that context, simulations of processes involving large strains, i.e. extrusion, can be performed with great robustness, using the same cloud of nodes during the whole simulation. Nevertheless, the NEM is computationally more expensive than FEM, as the shape functions can not be expressed analytically. Different natural neighbour interpolants (NEM shape functions) exist, i.e. Sibson shape functions, Laplace shape functions, or pseudo NEM shape functions. In this work, a comparison of these different interpolants in terms of CPU requirements and accuracy is proposed. Numerical investigations are performed through both linear 3D benchmark and complex extrusion simulations.

## 2 NATURAL NEIGHBOUR INTERPOLANTS

There exist various types of natural neighbour-based

interpolations, but they all are based on the construction of the Delaunay triangulation [1] of the cloud of points used to discretise the domain. The Delaunay triangulation of a set of points is the unique triangulation of the set that verifies the so-called empty circumcircle criterion. This means that no point of the set lies in the interior of a circle that passes through the three vertices of each triangle.

The dual structure of the Delaunay triangulation is the Voronoi diagram of the cloud [2]. For a given node $n_I$, its associated Voronoi cell is composed by all of the points which are closer to the node $n_I$ than to any other node. Formally,

$$T_I = \{x \in \mathbb{R}^n : d(x, x_i) < d(x, x_j) \forall j \neq i\} \quad (1)$$

where $T_I$ is the Voronoi cell and $d(\cdot, \cdot)$ represents the Euclidean distance. Two nodes whose Voronoi cells share one edge are called *natural neighbours* and hence the name of these interpolation schemes.

### 2.1 Sibson

The most extended form of natural neighbour-based interpolation is due to Sibson [3]. For the definition of Sibson interpolation it is necessary to previously introduce the concept of second order Voronoi cell $T_{II}$. It is defined as the locus of the points that have the node $n_I$ as the closest node and the node $n_J$ as the second closest node.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Areas of the second Voronoi cells</figcaption></figure>

If a new point $x$ (i.e. an integration point) is added to a given cloud of points, the Voronoi cells will be modified by its presence (see figure 1). Sibson [4] defined the natural neighbour coordinates of a point $x$ with respect to one of its neighbours $n_I$ as the ratio of the area (2D) or volume (3D) of the cell $T_{Ix}$ to the total area (2D) or volume (3D) of $T_x$:

$$ \phi_I^{\text{NEMS}}(x) = \frac{\text{Volume}(T_{Ix})}{\text{Volume}(T_x)} \qquad (2) $$

The resultant shape function obviously depends on the relative position of the nodes. An example of such shape function on a 2D regular mesh can be seen in figure 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Sibson interpolant</figcaption></figure>

## 2.2 Laplace

Laplace interpolants [5] are also based on the Voronoi diagram of the cloud of but involve computations on lower dimensional geometrical elements. Being $h_{xI}$ the distance between a node $n_I$ and a point $x$, and $m_{xI}$ the length (2D) or area (3D) of the common face of the cells defined by the node $n_I$ and the point $x$ (see figure 3), the Laplace function of node $n_I$ at the point $x$ is defined as:

$$ \phi_I^{\text{NEML}}(x) = \frac{\alpha_I(x)}{\alpha_T} \qquad (3) $$

where:

$$ \alpha_I(x) = \frac{m_{xI}}{h_{xI}} \qquad (4) $$

is called Laplace weight and

$$ \alpha_T = \sum_{I=1}^{n} \alpha_I \qquad (5) $$

is the sum of the weights of all $n$ neighbours nodes of $x$.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Parameters involved in Laplace interpolant computations</figcaption></figure>

Figure 4 shows the interpolant of a 2D regular mesh:

<figure><img src="image_5.png" /><figcaption>Fig. 4. Laplace interpolant</figcaption></figure>

## 3 PSEUDO NATURAL NEIGHBOUR INTERPOLANTS

In the context of Pseudo-NEM interpolants [6,7], the shape functions are computed using the Moving Least Square [8] framework, with particular weight functions, in order to match the properties of the NEM. For that purpose, some weight functions matching the support of the NEM shape functions support are computed as follows. The definition is given here in 2D, but is straightforward in 3D. Let a cone function which basis matches one of the Delaunay circle containing the point $x$, and where the projection of the node match the node $n_i$ (see figure 5).

<figure><img src="image_6.png" /><figcaption>Fig. 5. Conical function.</figcaption></figure>

The value of the conical function computed at point $x$ is given by:

$$ w(\mathbf{x}) = \frac{\|\mathbf{n}_i \mathbf{P}\| - \|\mathbf{n}_i \mathbf{x}\|}{\|\mathbf{n}_i \mathbf{P}\|}, \quad \mathbf{n}_i \mathbf{P} = -2 \left( \frac{\mathbf{c} \mathbf{n}_i \cdot \mathbf{n}_i \mathbf{x}}{\mathbf{n}_i \mathbf{x} \cdot \mathbf{n}_i \mathbf{x}} \right) \qquad (6) $$
In order to avoid overlapping of the cone functions whereas conserving the continuity of the weight function, a cone portion is associated with each Delaunay triangle connected to the node $n_i$ (see figure 6).

<figure><img src="image_7.png" /><figcaption>Fig. 6. Zone associated with a conical function.</figcaption></figure>

After introduction of (6) in the MLS framework, we obtain the shape functions depicted in figure 7 More details on the procedure can be found in [6,7].

<figure><img src="image_8.png" /><figcaption>Fig. 7. Pseudo NEM interpolant</figcaption></figure>

# 4 LINEAR ELASTIC COMPRESSION TEST

## 4.1 Problem definition

The first example consists in the compression of a cube made of an linear elastic material. A displacement of 0.01 on z direction is prescribed on the top face and the nodes of the bottom face are constrained in the z direction. Different clouds of nodes have been employed, both in regular and irregular distributions. In the FEM context, Delaunay triangulation has been used, and three integration points are used in each tetrahedron. All meshless simulations use the same integration points.

## 4.2 Results

### 4.2.a Computing time

Figure 8 shows the time needed to solve different problems. It is obvious that solving the problem using Sibson (NEMS) shape functions takes the highest times, increasing very fast. Furthermore, the

time needed for irregular distributions is higher than regular ones, since regular ones have less Delaunay tetrahedra and thus less integration points.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Computing time comparison</figcaption></figure>

Comparing the ratio of time needed to calculate the shape functions to the total time to solve the problem (figure 9), the pseudo NEM (PNEM) needs less than 10% of total time and decreases when the size of the problem increases. Laplace method (NEML) also decreases, but gets about 30% time. Again, Sibson is very expensive.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Shape calculation time/total time comparison</figcaption></figure>

# 5 EXTRUSION SIMULATION

## 5.1 Problem definition

The second problem consists in the simulation of a real extrusion problem. A regular mesh is generated at the beginning of the extrusion process. After some steps, the mesh becomes very distorted. The example involves 3021 nodes and 53781 integration points.

Mixed formulation with incompressible viscoplastic material model has been used. Only three iterations are needed for convergence.
## 5.2 Results

### 5.2.a Computing time

As can be seen in the table 2, only calculating Sibson shape functions have significance in the total time.

<table><caption>Table 2. shape function calculation time/total time for extrusion</caption><thead><tr><td></td><td>shape calc time</td><td>total time</td><td>proportion</td></tr></thead><tbody><tr><td>FEM</td><td>0.22</td><td>289.81</td><td>0.0759%</td></tr><tr><td>NEMS</td><td>17799.43</td><td>19611</td><td>90.7625%</td></tr><tr><td>NEML</td><td>70.94</td><td>1877.02</td><td>3.7794%</td></tr><tr><td>PNEM2</td><td>10.84</td><td>1814.12</td><td>0.5975%</td></tr></tbody></table>

### 5.2.b Accuracy

The following figures show the equivalent strain rate at a region close to the exit, where the mesh is highly distorted. The artificially values obtained with FEM due to the mesh distortion are not reproduced in the meshless methods. However, PNEM values are very low, probably due to the locking of the velocity field.

<figure><img src="image_11.png" /><figcaption>Fig. 10. Equivalent strain rate at the corner</figcaption></figure>

## 6 CONCLUSIONS

The calculation of the NEM and Pseudo NEM

interpolants takes more time than FEM shape functions. However, in complex calculation the extra time needed to calculate natural interpolants is relevant only when Sibson is used.

Due to the larger support of NEM shape functions, the matrices are not so empty, thus total calculation time increases. But when the high distortion of a mesh makes FEM unstable, meshless methods seem to be useful.

Despite Pseudo NEM involve less computational times than other NEM approximations, its efficiency seems quite poor compared to other reviewed shape functions. Nevertheless, it is worth noting, that the accuracy can be enhanced using special smoothing techniques [7], as well as special treatment of the boundary conditions. Another assumption is that the higher complexity of the shape functions probably necessitates a finer integration scheme. These improvements are currently being investigated.

## REFERENCES

1. B. Delaunay, Sur la Sphère Vide. A la mémoire de Georges Voronoi, In: Izvestia Akademii Nauk SSSR, Otdelenie Matematicheskii i Estestvennyka Nauk (1934) 7: 793-800.

2. G. M. Voronoi, Nouvelles Applications des Paramètres Continus à la Théorie des Formes Quadratiques. Deuxième Mémoire: Recherches sur les parallélloédres Primitifs, In: J. Reine Angew. Math. (1908) 134: 198-287.

3. R. Sibson, A Brief Description of Natural Neighbour Interpolation, In: Interpreting Multivariate Data, ed. V. Barnett, John Wiley (1981) 21-36.

4. R. Sibson, A Vector Identity for the Dirichlet Tesselation, In: Mathematical Proceedings of the Cambridge Philosophical Society (1980) 87: 151-155.

5. H. Hiyoshi, K. Sugihara, Two generalizations of an interpolant based on Voronoi diagrams, In: International Journal of Shape Modeling (1999) 5(2):219-231.

6. J. Yvonnet, F. Chinesta, A hybrid Moving Least Square natural element meshfree method for direct enforcement of boundary conditions and faster three dimensional computations, In: Third MIT conference on computational mechanics, Cambridge (2004).

7. J. Yvonnet, Nouvelles approaches sans maillage basées sur la méthode des éléments naturels pour la simulation des procédés de mise en forme, Phd Thesis (1994).

8. T. Belytschko, Lu G.G., Gu L., Element Free Galerkin methods. International Journal for Numerical Methods in Engineering (1994) 37: 229-256.