Arch. Comput. Meth. Engng.
Vol. 10, 4, 307-384 (2003)
Archives of Computational
Methods in Engineering
State of the art reviews
Overview and Recent Advances in Natural
Neighbour Galerkin M ethods
E. Cueto
Arag´on Institute for Engineering Research (I3A)
University of Zaragoza
Mar´ıa de Luna, 3. E-50018 Zaragoza. Spain
N. Sukumar
Department of Civil and Environmental Engineering
University of California, One Shield Avenue
Davis, CA 95616 USA
B. Calvo, J. Cego˜nino and M. Doblar´e
Arag´on Institute for Engineering Research (I3A)
University of Zaragoza
Mar´ıa de Luna, 3. E-50018 Zaragoza. Spain
mdoblare@posta.unizar.es
Summary
In this paper, a survey of the most relevant advances in natural neighbour Galerkin methods is presented. In
these methods (also known as natural element methods, NEM), the Sibson and the Laplace (non-Sibsonian)
interpolation schemes are used as trial and test functions in a Galerkin procedure. Natural neighbour-based
methods have certain unique features among the wide family of so-called meshless methods: a well-defined
and robust approximation with no user-defined parameters on non-uniform grids, and the ability to exactly
impose essential (Dirichlet) boundary conditions are particularly noteworthy.
A comprehensive review of the method is conducted, including a description of the Sibson and the Laplace
interpolants in two- and three-dimensions. Application of the NEM to linear and non-linear problems
in solid as well as fluid mechanics is studied. Other issues that are pertinent to the vast majority of
meshless methods, such as numerical quadrature, imposing essential boundary conditions, and the handling
of secondary variables are also addressed. The paper is concluded with some benchmark computations that
demonstrate the accuracy and the key advantages of this numerical method.
1INTRODUCTION
The finite element method (FEM), which was conceived in the 1950s, has now become the
most widely-used numerical method for computer simulation in academic research as well
as in engineering practice. Finite elements provide the possibility of handling arbitrary ge-
ometries with complex boundary conditions and nonlinear material behaviour with relative
ease—the rich mathematical framework to demonstrate convergence and adaptive refine-
ment and error estimation strategies have paved the way for its present success in modeling
and simulation. A practical drawback, however, is the need for regeneration of the mesh in
moving boundary and large deformation problems. This is often done by the analyst, and
is considered to be one of the most time-consuming tasks in a finite element analysis.
To overcome the difficulty associated with remeshing, the past decade has seen a
tremendous surge in the development of a family of Galerkin and collocation-based nu-
merical methods—these are known as particle, gridless, meshfree, or meshless methods.
For instance, smooth particle hydrodynamics [87], diffuse element method [91], element-
free Galerkin method (EFGM) [19], material point method (MPM) [112], reproducing ker-
nel particle method (RKPM) [79], finite point method [92], partition of unity method
[85, 11], h-p clouds [45], natural neighbour Galerkin methods, or natural element meth-
c
2003 by CIMNE, Barcelona (Spain). ISSN: 1134–3060 Received: November 2002
308 E. Cueto, N. Sukumar, B. Calvo, J. Cego˜nino and M. Doblar´e
ods (NEM) [24, 111], boundary node method (BNM) [89], meshless local Petrov-Galerkin
method (MLPG) [6], method of finite spheres [40], and finite cloud method [2] are some of
the most widely used techniques. We refer the interested reader to the review articles by
Belytschko and co-workers [17] and Li and Liu [77] for a detailed discussion and comparison
of different meshless and particle methods.
Meshless methods share the common characteristic of there being no explicit connec-
tivity information between nodes; the approximant is built in a process that is transparent
to the user. Most of these meshless methods emerged as a consequence of using exist-
ing interpolation or data approximation techniques, such as moving least squares (MLS)
approximation [73], as trial and test functions in a Galerkin procedure. Meshless approxi-
mants such as MLS do not interpolate nodal values and as a consequence, do not exactly
reproduce essential (Dirichlet) boundary conditions. This aspect as well as the topic of nu-
merical integration in meshless methods has attracted significant research interest in recent
years.
The natural element method (NEM) is a Galerkin-based meshless method that is built
upon the notion of natural neighbour interpolation. This interpolation scheme has very
striking properties, such as its strictly interpolating character, ability to exactly interpolate
piece-wise linear boundary conditions, and a well-defined and robust approximation with
no user-defined parameter on non-uniform grids.
In this paper, we review the most notable aspects of the NEM with emphasis on the
recent advances achieved by the authors in its application to solid as well as fluid mechan-
ics. The outline of the paper is as follows: in Section 2, some well-known computational
geometric constructs, and a natural neighbour-based (Sibson and Laplace) interpolation
schemes are described in detail. The imposition of essential boundary conditions is anal-
ysed in Section 3, and issues concerning errors in numerical integration and the description
of a stabilized nodal quadrature scheme are presented in Section 4. A mixed approximation
with application to incompressible and nearly incompressible elasticity is discussed in Sec-
tion 5. In Section 6 the ability of the method to handle fluid mechanics problems with free
surfaces (e.g., mould filling) is addressed. In Section 7, the Laplace interpolant is used to
construct a finite difference scheme on unstructured grids. Applications in fluid mechanics
and biomechanics are pursued in Section 8, and finally in Section 9, the main conclusions
that are drawn from this study are mentioned.
2 NATURAL NEIGHBOUR-BASED INTERPOLANTS
Sibson [103, 104] introduced the notion of natural neighbour interpolation for data ap-
proximation and smoothing, and further investigations on the construction and properties
of the natural neighbour (or Sibson) interpolant were carried out by Farin [50] and Piper
[95]. Recently, renewed interest has sparked in the computational geometry community on
the Sibson interpolant and generalization of Voronoi-based interpolation schemes that are
based on natural neighbours [59, 60, 61, 20].
2.1 Voronoi Diagram and Delaunay Triangulation
Consider a bounded domain Ω in d-dimensions described by a set N of M scattered nodes:
N = {n
1
,n
2
,... ,n
M
}. The Voronoi diagram [9] V(N)ofthesetN is a sub-division of
the domain into regions V (n
I
), such that any point in V (n
I
)isclosertonoden
I
than to
any other node n
J
∈ N (J = I). The region V (n
I
) (first-order Voronoi cell) for a node n
I
within the convex hull is a convex polygon (polyhedron) in
R
2
(
R
3
):
V (n
I
)={x ∈
R
d
: d(x, x
I
) <d(x, x
J
) ∀J = I}, (1)
Overview and Recent Advances in Natural Neighbour Galerkin Methods 309
where d(x
I
, x
J
) is an appropriate distance function (usually the standard Euclidean distance
is used) between x
I
and x
J
.
The dual of the Voronoi diagram, the Delaunay tessellation, is constructed by connecting
nodes that have a common (d–1)-dimensional Voronoi facet. Given any nodal set N,the
Voronoi diagram is unique, whereas the Delaunay tessellation is not—a simple example is
the triangulation of a square where choosing either diagonal leads to two valid Delaunay
triangulations. In Figure 1a, the Voronoi diagram and the Delaunay triangulation are
shown for a nodal set consisting of seven nodes (M = 7). A Voronoi vertex and an edge
are also indicated in Figure 1a. An important property of Delaunay triangles is the empty
circumcircle criterion [75]—if DT (n
J
,n
K
,n
L
) is any Delaunay triangle of the nodal set N,
then the circumcircle of DT contains no other nodes of N. In Figure 1b, the Delaunay
circumcircles for three triangles are shown.
Consider now the introduction of a point p with coordinates x ∈
R
2
into the domain
Ω (Figure 1b). The Voronoi diagram V(n
1
,n
2
,... ,n
M
,p) or equivalently the Delaunay
triangulation DT (n
1
,n
2
,... ,n
M
,p)fortheM nodes and the point p is constructed. Now,
if the Voronoi cell for p and n
I
have a common facet (segment in
R
2
and polygon in
R
3
),
then the node n
I
is said to be a natural neighbour of the point p [103]. The Voronoi cells
for the point p and its natural neighbours are shown in Figure 1c, together with the convex
hull of the set of points.
2.2 Sibson Interpolant
The natural neighbour (Sibson) interpolant was introduced by Sibson [103]. For simplicity
and ease of exposition, we restrict our attention to 2-dimensions, altough every concept is
easily extended to d>2. The definition of the Voronoi diagram (1st-order) appears in
Eq. (1). By a similar extension, one can construct higher order (k-order, k>1) Voronoi
diagrams in the plane. Of particular interest in the present context is the case k =2,
which is the second-order Voronoi diagram. The second-order Voronoi diagram of the set
of nodes N is a sub-division of the plane into cells V
IJ
, such that V
IJ
is the locus of all
points that have n
I
as the nearest neighbour, and n
J
as the second nearest neighbour. The
second-order Voronoi cell V
IJ
(I = J) is defined as [103]
V
IJ
= {x ∈
R
2
: d(x, x
I
) <d(x, x
J
) <d(x, x
K
) ∀K = I, J}. (2)
In order to quantify the neighbours for a point p with coordinate x =(x
1
,x
2
)thatis
inserted into the tessellation, Sibson [103] used the concept of second-order Voronoi cells,
and thereby introduced natural neighbours and natural neighbour coordinates. Natural
neighbour coordinates (shape functions) are used as the interpolating functions in natural
neighbour (Sibson) interpolation, and as trial and test functions in a Galerkin implemen-
tation for the solution of partial differential equations (PDEs) [106]. Consider Figure 2a,
where a point p is inserting into a tessellation. The natural neighbour shape function of p
with respect to a natural neighbour I is defined as the ratio of the area of the second-order
Voronoi cell (A
I
) to the total area of the first-order Voronoi cell of p (A):
φ
I
(x)=
A
I
(x)
A(x)
,A(x)=
n
J=1
A
J
(x), (3)
where n =5forthepointp in Figure 2a. In 3-d, the Sibson shape function is defined as
the ratio of polyhedral volumes.
The derivatives of the Sibson shape functions are obtained by differentiating Eq. (3):
φ
I,j
(x)=
A
I,j
(x) − φ
I
(x)A
,j
(x)
A(x)
(j =1, 2), (4)
310 E. Cueto, N. Sukumar, B. Calvo, J. Cego˜nino and M. Doblar´e
P
i
Delaunay triangle
Voronoi edge
Voronoi cell
Voronoi
vertex
(a)
1
p
2
3
4
5
6
7
(b)
1
p
2
3
4
5
6
Convex hull
7
(c)
Figure 1. (a) Voronoi diagram and Delaunay triangulation; (b) Delaunay circum-
circles and (c) Natural neighbours (filled circles) of inserted point p
where Eq. (3) has been used to arrive at the above expression. If the point x → x
I
,then
φ
I
(x) = 1 and all other shape functions are zero. Therefore, the properties of positivity,
interpolation, and partition of unity are straightforward [110]:
0 ≤ φ
I
≤ 1,φ
I
(x
J
)=δ
IJ
,
n
I=1
φ
I
(x)=1. (5)
Overview and Recent Advances in Natural Neighbour Galerkin Methods 311
1
p
2
3
6
7
A(p)
A
1
(a)
1
p
2
3
6
7
s
2
s
3
s
6
s
7
s
1
h
2
h
3
h
1
h
6
h
7
(b)
Figure 2. Natural neighbour-based interpolants. (a) Sibson interpolant; and (b)
Laplace interpolant
Natural neighbour shape functions also satisfy the local coordinate property [103], namely
x =
n
I=1
φ
I
(x)x
I
, (6)
and hence in conjunction with the partition of unity property, it is readily derived that the
Sibson interpolant can exactly represent any linear field which is known as linear complete-
ness in the finite element literature [18].
2.3 Laplace Interpolant
We trace the roots of the Laplace interpolant with a view to summarize and unify some
of the previous developments. As first noted in [107], Belikov and co-workers [14, 15, 16]
as well as Hiyoshi and Sugihara [105, 59, 60] independently proposed a natural neighbour-
based interpolant that was different from the Sibson interpolant. The former scientists who
are in the field of data approximation and PDEs referred to the new interpolant as the non-
Sibsonian interpolant, whereas Hiyoshi and Sugihara (computational geometers) coined it
as the Laplace interpolant. Both groups recognized its connection to the Laplace equation,
and used very similar approaches in delineating many of its properties. In this paper, we
choose the name of Laplace interpolant. In [14, 15, 16], tools from vector calculus and data
approximation theory are used to investigate the Laplace interpolant; Hiyoshi and Sugihara
[105, 60] use the Minkowski theorem for convex polytopes [56] (Gauss’s divergence theorem)
and geometric-based theorems to study the Sibson and Laplace interpolants. They view
the Laplace and Sibson interpolants as particular instances of Voronoi-based interpolants,
and present a framework for generating a hierarchy of natural-neighbour based interpolants
with increasing order of continuity at non-nodal locations. The Sibson interpolant can
be obtained by the Voronoi-based integration of the Laplace interpolant and therein lies
a reason why the Sibson interpolant is smoother than the Laplace interpolant [60]. By
