# Radial Basis Point Interpolation Method with Reordering Gauss Domains for 2D Plane Problems

Authors: Shi-Chao Yi; Fu-jun Chen; Lin-Quan Yao
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We present novel Gauss integration schemes with radial basis point interpolation method (RPIM). These techniques define new Gauss integration scheme, researching Gauss points (RGD), and reconstructing Gauss domain (RGD), respectively. The developments lead to a curtailment of the elapsed CPU time without loss of the accuracy. Numerical results show that the schemes reduce the computational time to 25% or less in general.


---

## Body

## 1. Introduction

In the past two decades, the development and application of meshfree methods have attracted much attention. One of the reasons is the versatility of meshfree methods for complex geometry of solids and flexibility for different engineering problems [1]. Element-free Galerkin (EFG) method, which is originated by Belytschko et al. [2–5], is one of the most widely used meshfree methods. The key advantage of EFG method is that only nodal data is required and no element connectivity is needed, when moving least-squares (MLS) interpolation is used to construct trial and test functions. It is currently widely used in computational mechanics and other areas, such as by Sarboland and Aminataei [6] for nonlinear nonhomogeneous Burgers Equation and Pirali et al. [7] for crack discontinuities problem. In the meantime, some new techniques are used to improve the performance of the MLS, complex variable moving least-squares [8–10] and improved complex variable moving least-squares [11], the moving least-squares with singular weight function [12], and so forth. But shape functions constructed by MLS interpolation do not possess Kronecker delta function property; the treatment of essential boundary conditions is one of the typical drawbacks. Thus, many special techniques have been proposed to impose essential boundary conditions [13–15], such as point collocation [13], Lagrange multipliers [2], singular weighting functions [14], and penalty method [15]. None of these methods is fully satisfactory, as they still need additional efforts to enforce essential boundary conditions.In order to totally eliminate the drawback associated with EFG method for imposing essential boundary conditions, Liu and Gu have developed the point interpolation methods (PIM) by using polynomial basis or/and radial basis function (RBF) [16–19]. However, singularity may occur if arrangement of the nodes is not consistent with the order of polynomial basis, while the inverse of moment matrix always exists for arbitrarily scattered nodes with radial basis [20]. In this paper, we mainly care about the radial point interpolation method (RPIM) studied by Wang and Liu in [21], which is based on the global weak form.On the other hand, high computational cost is still one of the main drawbacks in meshfree method with Galerkin weak form. In the RPIM method for a given computational point, ann
×
n linear system (the coefficient matrix is called the moment matrix) should be solved to construct shape functions if radial basis functions are selected, where n is the number of nodes in the support domain. Furthermore, if m monomial basis functions are added, the linear system will be extended to (
n
+
m
)
×
(
n
+
m
). This is very time-consuming especially for the meshless methods based on weak forms, where integration should be employed and the number of computational points is much bigger than the number of nodes.In fact, common meshless method is based on the numerical integration for Gauss domains with RPIM. It is a tedious process when the integration points are extremely more than the nodes. In practical work, some Gauss point has the same nodes as some neighboring computational points. Thus, they only need one shape function with them. A process will be needed to find these nodes. This scheme is called researching Gauss points method (RGP). But the storage of searching program will increase quickly with the increasing number of nodes and Gauss points. Thus, more computational time is spent. To avoid this, reconstructing Gauss domains (RGD) is presented. All computational points share one same influence domain (the same nodes) in the domain. This scheme offsets the RGP’s disadvantage without the searching program. The reason why we call it researching Gauss points method is that the RGP method needs a searching nodes program same as in the RPIM. The two techniques are collectively named reordering Gauss domain (RG) method.The remainder of the paper is arranged as follows. In Section2, the radial basis point interpolation method (RPIM) is presented to get the shape functions. In Section 3, a Galerkin weak form and its numerical algorithm are studied for 2D solid mechanics problems. The reordering Gauss domains methods are then constructed in Section 4. In Section 5, several examples are presented to show the effectiveness of the proposed method and some parameters’ performances of the proposed method are also investigated. Finally, we end this paper with some conclusions in Section 6.


## 2. Radial Basis Point Interpolation Method (RPIM)

The RPIM interpolationu
h
(
x
) of u, for all x
∈
Ω
x, can be defined by
(1)
u
h
(
x
)
=
∑
i
=
1
n
‍
R
i
(
x
)
a
i
+
∑
j
=
1
m
‍
p
j
(
x
)
b
j
with the constraint condition
(2)
∑
i
=
1
n
‍
p
j
(
x
i
)
a
i
=
0
,
j
=
1,2
,
…
,
m
,
where R
i
(
x
) is the radial basis function (RBF), n is the number of nodes in the neighborhood of x, p
j
(
x
) is the monomial in the space coordinates p
T
(
x
)
=
[
1
,
x
,
y
], m is the number of polynomial basis functions, and coefficients a
i and b
j are interpolation constants. The variable x of the radial basis function R
i
(
x
) is the distance between the interpolation point x and a node x
i. It is necessary to construct the interpolation function here to solve the equations.There are a number of types of radial basis functions. Characteristics of radial basis functions have been widely investigated in [20, 22]. In this paper, the following multiquadrics (MQ) radial basis function is used:
(3)
R
i
(
x
)
=
r
i
2
+
C
2
q
,
where r
i
=
x
-
x
i
2
+
y
-
y
i
2 and C and q are two parameters. C is defined as
(4)
C
=
α
0
d
i
,
where α
0 is a dimensionless coefficient and d
i is a parameter of the nodal distance. For regularly distributed nodal case, d
i is the shortest distance between the node x
i and its neighbor nodes. Effects of α
0 and q have been studied in detail in [22]. In static analysis of 2D solid problems, it has been found that α
0
=
1.0 and q
=
1.03 lead to good results. Hence, these numbers are used in this paper.The second term in (1) consists of polynomials. To ensure invertible interpolation matrix of RBF, the polynomial added into the RBF cannot be arbitrary. A low-degree polynomial is often needed to augment RBF to guarantee the nonsingularity of the moment matrix. In addition, the linear polynomial added into the RBF can also ensure linear consistency and improve the interpolation accuracy [22]. Thus, the linear polynomial is added into the MQ RBF in the following discussion.Coefficientsa
i and b
j in (1) can be determined by enforcing that (1) and (2) be satisfied at the n nodes surrounding point x. Equations (1) and (2) can be rewritten in the matrix form:
(5)
G
a
0
≜
R
0
P
0
P
0
T
0
a
b
=
u
s
0
,
where
(6)
a
0
T
=
a
T
b
T
T
=
a
1
,
a
2
,
…
,
a
n
,
b
1
,
b
2
,
b
3
,
u
s
T
=
u
1
,
u
2
,
…
,
u
n
,
R
0
=
R
1
(
r
1
)
R
2
(
r
1
)
⋯
R
n
(
r
1
)
R
1
(
r
2
)
R
2
(
r
2
)
⋯
R
n
(
r
2
)
⋮
⋮
⋱
⋮
R
1
(
r
n
)
R
2
(
r
n
)
⋯
R
n
(
r
n
)
n
×
n
,
P
0
T
=
1
1
⋯
1
x
1
x
2
⋯
x
n
y
1
y
2
⋯
y
n
3
×
n
.
The matrix R
0 is symmetric, and thus the matrix G is symmetric. Then, the interpolation equation (1) is finally expressed as
(7)
u
h
(
x
)
=
R
T
x
p
T
x
G
-
1
u
s
=
Φ
R
(
x
)
u
s
,
where the shape function Φ
(
x
) is defined by
(8)
Φ
(
x
)
=
ϕ
1
x
,
ϕ
2
x
,
…
,
ϕ
n
x
=
R
T
x
p
T
x
G
-
1
.
And R
T
(
x
)
=
[
R
1
(
x
)
,
R
2
(
x
)
,
…
,
R
n
(
x
)
], p
T
(
x
)
=
[
1
,
x
,
y
]. It can be found from the above discussion that RPIM passes through the nodal values. Therefore, RPIM shape functions given in (8) satisfy the Kronecker delta condition. Thus,
(9)
ϕ
i
(
x
j
)
=
δ
i
j
=
1
,
i
=
j
,
0
,
i
≠
j
.


## 3. Variational Form of 2D Plane Problem

Consider the 2D problem of the deformation of a linear elastic medium from an undeformed domainΩ, enclosed by Γ:
(10)
∇
·
σ
+
b
=
0
,
in
Ω
,
where σ is the stress tensor corresponding to the displacement field u and b is a body force vector, and boundary conditions are as follows:
(11)
σ
·
n
=
t
¯
on
⁡
Γ
t
,
u
=
u
¯
on
⁡
Γ
u
,
in which t
¯ and u
¯ are prescribed tractions and displacements, respectively, on the traction boundary Γ
t and on the displacement boundary Γ
u and n is the unit outward normal matrix to the boundary Γ
t.Using the standard principle of minimum potential energy for (10)-(11), that is, to find u
∈
H
1
ω
3 such that
(12)
Π
=
1
2
∫
Ω
‍
ε
T
D
ε
d
Ω
-
∫
Ω
‍
u
T
·
b
d
Ω
-
∫
Γ
t
‍
u
T
·
t
¯
d
Γ
is stationary, where H
i
(
ω
) denotes the Sobolev space of order i, ε and σ
=
D
ε are strain-stress vectors, and D is the strain-stress matrix. RPIM equation (1) is used to approximate the displacements in the Galerkin procedure. Then we can get
(13)
u
=
u
x
u
y
=
∑
i
=
1
k
‍
ϕ
i
0
0
ϕ
i
u
x
i
u
y
i
=
∑
i
=
1
k
‍
Φ
i
u
i
.
Substituting (13) into (12) leads to the following total potential energy in the matrix form:
(14)
Π
=
1
2
u
T
·
K
·
u
-
u
T
·
f
,
and invoking δ
Π
=
0 results in the following linear system of u:
(15)
K
u
=
f
,
in which the stiffness matrix K is built from 2
×
2 matrices K
i
j and the right-hand side vector f is built from the 2
×
1 vectors f
i. These matrices and vectors are defined by
(16)
K
i
j
=
∫
Ω
‍
B
i
T
D
B
j
d
Ω
,
i
,
j
=
1,2
,
…
,
k
,
f
i
=
∫
Ω
‍
ϕ
i
b
d
Ω
+
∫
Γ
t
‍
ϕ
i
t
¯
d
Γ
,
i
=
1,2
,
…
,
k
,
where
(17)
B
i
=
ϕ
i
,
x
0
0
ϕ
i
,
y
ϕ
i
,
y
ϕ
i
,
x
,
D
=
E
1
-
ν
2
1
ν
0
ν
1
0
0
0
1
-
ν
2
(
for
plane
stress
problem
)
.Whether the RPIM method or other meshless methods based on global weak form, background cells are necessary to obtain the numerical integration of (16). Different from the finite element method (FEM), which uses the same nodes for both interpolation and numerical integration, background cells in meshless methods are independent of the interpolations. In this paper, quadrilateral cells and Gauss quadrature are used for the numerical integration.


## 4. Reconstructing Gauss Domains Methods

A weak form, in contrast to a strong form (collocation method in general), requires weaker consistency on the assumed field functions. The consistency requirement on the assumed functions for field variables is very different from the strong form. For a2
kth-order differential governing system equation, the strong formulation requires a consistency of the 2
kth order, while the weak formulation requires a consistency of only the kth order. But the meshfree method usually uses the integral representation of field variable functions for solving strong form system equations and the numerical integration is extremely time-consuming. Numerical accuracy mainly depends on the number of Gauss points in the corresponding domain; the more the Gauss points, the better the results in general. Thus, one of the most time-consuming steps in the meshless method is the construction of shape functions, since for every point of interest a linear system should be computed. As discussed in Section 2, using radial basis functions to construct shape functions, an n
×
n linear system should be computed for every computational point. If m monomials are added, an (
n
+
m
)
×
(
n
+
m
) linear system should be solved. This is very time-consuming especially for meshless methods based on weak form, where a large number of integration points are used. In this section, we propose the researching Gauss points (RGP) and reconstructing Gauss domains (RGD) methods, which are together named reordering Gauss domains (RG) methods and can partly reduce the computation cost for the meshless methods compared with the RPIM.


### 4.1. Researching Gauss Points (RGP) with the Same Nodes

In the RPIM approximation, shape function consists of two parts:R
T
(
x
)
p
T
(
x
) and G
-
1   (8). The time consumption of R
T
x
p
T
x is less than that of G
-
1, because the computational complexities are o
(
n
+
m
) and o
(
(
n
+
m
)
3
), respectively.Every Gauss point has its ownR
T
(
x
)
p
T
(
x
) and is different from the rest, but it may have the same G
-
1 as the other. Thus, we need a storage place, containing the public nodes and the Gauss points’ information. The data is obtained by searching the Gauss points. Thus, the method is called researching Gauss points (RGP) method.The red area represents the nodes set (involved red nodes∘) and the blue area represents the Gauss points set (involved blue gauss points *) in Figure 1(a). It should be pointed out that one-to-one mapping exists between the Gauss points sets and the nodes points sets (can be verified by logical deduction). Thus, the key work is how to get the relationship. First of all, we construct N
n
*
N
g matrix S
=
[
s
i
j
], where N
n
,
N
g denote the number of the nodes and the number of the Gauss points, respectively. The element S
=
[
s
i
j
] is δ function:
(18)
s
i
j
=
1
the
i
th
nodes
locate
in
the
support
domain
of
the
j
th
gauss
point
;
0
otherwise
.The distribution of the nodes and the Gauss points: find the 1-1 mapping (a); assign the 1-1 mapping (b).
(a)
(b)Thus,S is written as
(19)
S
=
⋯
⋯
⋯
⋯
⋯
1
1
1
⋯
0
⋯
⋯
⋯
⋯
⋯
1
1
1
⋯
0
⋯
⋯
⋯
⋯
⋯
N
n
*
N
g
.If thejth row elements are the same as the kth row elements, we confirm the jth and kth Gauss points have the same basis functions R
(
x
j
)
R
(
x
k
). It is a rule, which gives one-to-one mappings between the two types of sets by comparing the rows of S. Some codes (isequal.m, unqiue.m, etc.) could be used to find these mappings in MATLAB. The process of one-to-one mapping is also called the searching program distinguishing the searching nodes program in advance. For the sake of simplicity, the RGP method is divided into two parts: Part I (the researching program) and the rest.


### 4.2. Reconstructing Gauss Domains (RGD)

The searching program saves the time consumption in Gauss integration, but it needs extra time and storage for the searching process. With the rapid increase of the nodes and the Gauss points, the searching process costs also grow exponentially. Then, we present the reconstructing Gauss domains (RGD) method without the searching process. In the fixed Gauss domain, all Gauss points are mandatorily contacted with some neighbouring nodes. These nodes lie in the support domain of the Gauss domain’s centroid. That is, every Gauss domain has its own basis functions orG
-
1 corresponding to the centroid.A center point (×) of the Gauss domain is added in Figure 1(b) compared with Figure 1(a). Without the searching process, we mandatorily consider that the center point is seen by not only the barycentre of the Gauss domain, but also the barycentre of the domain involving some nodes. These nodes locate in the support domain with α on the barycentre point. Thus, the one-to-one mappings are ascertained by the corresponding barycentre points.


## 4.1. Researching Gauss Points (RGP) with the Same Nodes

In the RPIM approximation, shape function consists of two parts:R
T
(
x
)
p
T
(
x
) and G
-
1   (8). The time consumption of R
T
x
p
T
x is less than that of G
-
1, because the computational complexities are o
(
n
+
m
) and o
(
(
n
+
m
)
3
), respectively.Every Gauss point has its ownR
T
(
x
)
p
T
(
x
) and is different from the rest, but it may have the same G
-
1 as the other. Thus, we need a storage place, containing the public nodes and the Gauss points’ information. The data is obtained by searching the Gauss points. Thus, the method is called researching Gauss points (RGP) method.The red area represents the nodes set (involved red nodes∘) and the blue area represents the Gauss points set (involved blue gauss points *) in Figure 1(a). It should be pointed out that one-to-one mapping exists between the Gauss points sets and the nodes points sets (can be verified by logical deduction). Thus, the key work is how to get the relationship. First of all, we construct N
n
*
N
g matrix S
=
[
s
i
j
], where N
n
,
N
g denote the number of the nodes and the number of the Gauss points, respectively. The element S
=
[
s
i
j
] is δ function:
(18)
s
i
j
=
1
the
i
th
nodes
locate
in
the
support
domain
of
the
j
th
gauss
point
;
0
otherwise
.The distribution of the nodes and the Gauss points: find the 1-1 mapping (a); assign the 1-1 mapping (b).
(a)
(b)Thus,S is written as
(19)
S
=
⋯
⋯
⋯
⋯
⋯
1
1
1
⋯
0
⋯
⋯
⋯
⋯
⋯
1
1
1
⋯
0
⋯
⋯
⋯
⋯
⋯
N
n
*
N
g
.If thejth row elements are the same as the kth row elements, we confirm the jth and kth Gauss points have the same basis functions R
(
x
j
)
R
(
x
k
). It is a rule, which gives one-to-one mappings between the two types of sets by comparing the rows of S. Some codes (isequal.m, unqiue.m, etc.) could be used to find these mappings in MATLAB. The process of one-to-one mapping is also called the searching program distinguishing the searching nodes program in advance. For the sake of simplicity, the RGP method is divided into two parts: Part I (the researching program) and the rest.


## 4.2. Reconstructing Gauss Domains (RGD)

The searching program saves the time consumption in Gauss integration, but it needs extra time and storage for the searching process. With the rapid increase of the nodes and the Gauss points, the searching process costs also grow exponentially. Then, we present the reconstructing Gauss domains (RGD) method without the searching process. In the fixed Gauss domain, all Gauss points are mandatorily contacted with some neighbouring nodes. These nodes lie in the support domain of the Gauss domain’s centroid. That is, every Gauss domain has its own basis functions orG
-
1 corresponding to the centroid.A center point (×) of the Gauss domain is added in Figure 1(b) compared with Figure 1(a). Without the searching process, we mandatorily consider that the center point is seen by not only the barycentre of the Gauss domain, but also the barycentre of the domain involving some nodes. These nodes locate in the support domain with α on the barycentre point. Thus, the one-to-one mappings are ascertained by the corresponding barycentre points.


## 5. Numerical Experiments

In this section, several numerical examples are selected to demonstrate the applicability of the RG meshless methods. The numerical results for these examples are compared with the analytical solutions and the RPIM solution. Square support domains are used for calculations in the present paper.d
c stands for the average node distance. The MQ radial basis function and the linear polynomial are used to construct shape functions. All runs are performed in MATLAB 7.0 on an Intel Pentium 4 (2 GB RAM) Windows XP system.


### 5.1. Patch Test

The first numerical example is the standard patch test shown in Figure2. The patch test consists of thirteen nodes including five interior nodes. A 2
×
2 rectangular background mesh is used for numerical integration and 4
×
4 Gauss points are used in each mesh. In this patch test, the displacements are prescribed on all outside boundaries by a linear function of x and y on the patches of the dimension. The parameters are taken as E
=
1.0 and ν
=
0.3. The linear displacement functions are u
x
=
0.6
x and u
y
=
0.6
y. RG methods of the patch test require that the displacement of any interior node be given by the same linear function and that the strains and stresses be constant in the patch. The RG methods pass the patch test when linear polynomials are added (m
=
3). However, when polynomial term is not included (m
=
0), the patch test does not easily pass. The same is for the RPIM method; see [21]. Computed results at interior nodes for the RG methods and the RPIM method and the exact results are listed in Table 1.Table 1
Results at interior nodes located irregularly in the 13 nodes’ patch.
Internal node
Coordinates
Exact
RGP/RPIM
RGD
5
(1.0, 1.0)
(
u
x
,
u
y
)
(0.600, 0.600)
(0.601, 0.599)
(0.601, 0.599)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.853, 0.853, 0)
(0.853, 0.853, 0)
10
(0.65, 1.0)
(
u
x
,
u
y
)
(0.390, 0.600)
(0.391, 0.599)
(0.391, 0.599)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.861, 0.853, 0.002)
(0.861, 0.853, 0.002)
11
(0.70, 1.5)
(
u
x
,
u
y
)
(0.420, 0.900)
(0.422, 0.897)
(0.421, 0.897)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.856, 0.855, 0.001)
(0.856, 0.855, 0.001)
12
(1.3, 1.2)
(
u
x
,
u
y
)
(0.780, 0.720)
(0.779, 0.718)
(0.779, 0.718)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.846, 0.852, −0.002)
(0.846, 0.852, −0.002)
13
(1.2, 0.6)
(
u
x
,
u
y
)
(0.720, 0.360)
(0.720, 0.359)
(0.720, 0.359)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.854, 0.856, 0)
(0.854, 0.856, 0)Figure 2
A patch with 13 nodes.


### 5.2. Cantilevered Beam

The second example is a cantilever beam problem; see Figure3. Consider a beam of length L and height H subjected to traction at the free end. The beam has a unit thickness, and thus a plane stress problem is considered here. The closed-form solution is available for parabolic traction of force P:
(20)
u
x
=
P
y
6
E
I
6
L
-
3
x
x
+
2
+
ν
y
2
-
H
2
4
,
u
y
=
-
P
6
E
I
3
ν
y
2
L
-
x
+
4
+
5
ν
H
2
x
4
+
3
L
-
x
x
2
,
where the moment of inertia I of the beam is given by I
=
H
3
/
12. The stresses corresponding to the above displacements are
(21)
σ
x
=
P
L
-
x
y
I
,
τ
x
y
=
-
P
2
I
H
2
4
-
y
2
,
σ
y
=
0
.
The beam parameters are taken as E
=
3.0
×
1
0
7,  ν
=
0.3,  H
=
12,  L
=
48, and P
=
1000 in computation.Figure 3
Cantilever beam.To evaluate effects of various parameters, we use the two-norm relative errorsr
u and r
σ for displacement and stress, respectively. They are defined by
(22)
r
u
=
u
-
u
^
2
u
2
,
r
σ
=
σ
^
-
σ
2
σ
2
,
where u
^ and u are the approximate and exact solution of displacements and σ
^ and σ are the approximate and exact value of stresses.


#### 5.2.1. Effect of Irregular Node Distribution

Node distributions with 325 irregular nodes and 481 irregular nodes are shown in Figure4. A background mesh is necessary to obtain the numerical integrations of (15). This mesh is independent of nodes for interpolations, while FEM uses the same nodes for both interpolation and numerical integration. For this problem, 10
×
6 and 15
×
6 background cells are used for the 325 irregular nodes’ problem and 481 irregular nodes’ problem, respectively. In each cell, a 4
×
4 Gauss quadrature is used to evaluate the stiffness matrix. Figure 5(a) shows a comparison of the analytical solution and numerical solution of u
y along y
=
0 for this problem with 325 irregular nodes. The plot shows an excellent agreement between the analytical and numerical results for each method when an irregular node distribution is used. Figure 5(b) is the sectional distribution of shear stress τ
x
y along the x
=
24 section for this problem with 325 irregular nodes. The closed-form solution is also plotted for comparison.Meshless models with irregular data points: 325 nodes (a), 481 nodes (b).
(a)
(b)Results for 325 irregular nodes: (a) displacement and (b) shear stress.
(a)
(b)Computational results including error and computational cost are listed in Table2 for the cantilevered beam problem with 325 irregular nodes and 481 irregular nodes. From Figure 5 and Table 2, we can find that the RG methods have better accuracy and are less time-consuming than the RPIM method.Table 2
Computational results for irregular node distribution.
Irregular nodes
RPIM
RGP(Part I + the rest)
RGD
325
r
u
0.0046
0.0013
r
σ
0.0793
0.0582
CPU
3.034
2.917 = 1.424 + 1.493
1.162
481
r
u
0.0080
0.0019
r
σ
0.0824
0.0472
CPU
4.521
6.468 = 4.354 + 2.114
1.668


#### 5.2.2. Effect of the Size of Support Domain of RGD

As discussed in Section4.2, the centroid of every Gauss cell has a support domain. The range of the support domain α is an important parameter. The 325 regular nodes’ distribution (25
×
13: 25 nodes in the x direction and 13 nodes in the y direction) is used to study the effect of the size of the support domain for the RG methods.In Table3, we list the numerical results for the RGD method with different α. To compare with the RPIM method, we also list the computational results for the RPIM method with different radiuses of the support domain. From Table 3, we can see that the RGD method has better steadiness, computational efficiency, and computational accuracy than the RPIM method. Particularly, with α
=
4.5, the elapsed CPU time is reduced to 12.7
% (see the data in the boxes) of the previous one. This shows that the size of the support domain has increased and its CPU time takes up a larger proportion of the entire time, in which the consumption surged from 48
% to 76
% of the entire cost. The searching process has become the major expenditure of the entire program. From Table 3, we can find that the RGD method has a stable convergency with the increase of the α, but it does not appear in the RPIM. In this paper, we take the suitable α to be 4.Table 3
Effect of the size of support domain of RGD (Example5.2).
α
RGD
RPIM and RGP
r
u
r
σ
CPU
r
u
r
σ
RPIM CPU
RGP CPU (Part I + the rest)
2.5
0.0513
0.2033
0.413
0.0074
0.0522
1.808
2.151
=
1.023
+
1.128
3.0
0.0129
0.1000
0.646
0.0031
0.0620
3.003
2.537
=
1.131
+
1.406
3.5
0.0109
0.0945
0.733
0.0031
0.0305
4.752
3.277
=
1.798
+
1.479
4.0
0.0006
0.0285
1.182
0.0002
0.0533
7.135
4.071
=
2.608
+
1.463
4.5
0.0004
0.0274
1.336
0.0021
0.0191
10.303
5.195
=
3.720
+
1.475
5.0
0.0002
0.0202
2.112
0.0004
0.0294
14.599
6.451
=
4.911
+
1.540


#### 5.2.3. Convergence Study

For convergence studies, five different regular node distributions with 52 (13
×
4: 13 nodes in the x direction and 4 nodes in the y direction) nodes, 175 (25
×
7) nodes, 370 (37
×
10) nodes, 637 (49
×
13) nodes, and 976 (61
×
16) nodes and five different irregular node distributions are considered. The size of the support domain is taken as α
=
4. The computational results are listed in Tables 4 and 5. In Tables 4 and 5, we also list the computational results of the RPIM method. The convergence curves with different node distributions are plotted in Figure 5. In Figure 5, h is the maximum size of node arrangement. From Tables 4 and 5 and Figure 5 we can see that the present method possesses good convergence and less computational times than the RPIM method in general. However, the searching process consumption overtakes 88
% of the entire process in the 976 regular/irregular nodes, which become cumbersome. The oscillatory behavior of the convergence curve of RGD method exists in regular node, and oscillatory behaviors appear in the two node distributions of the RPIM method. Figure 5 shows that RGD method is steadier than the RPIM method. The reason of oscillatory appearance is that the accuracy of most meshless methods is closely related to the integration scheme, the number of Gauss points, the radiuses of the support domain and influence domain, and so on. Thus, it is a difficult task to get the best accuracy for all cases.Table 4
Convergence study on regular nodes for Example5.2.
Regular nodes
13
×
4
25
×
7
37
×
10
49
×
13
61
×
16
RGD
r
u
0.0384
0.0063
0.0035
0.0025
0.0015
r
σ
0.1891
0.0399
0.0676
0.0668
0.0182
CPU
0.440
0.839
1.972
2.730
5.420
RGP/RPIM
r
u
0.0358
0.0028
0.0024
0.0196
0.0240
r
σ
0.2057
0.0522
0.0559
0.1434
0.1738
RPIM CPU
1.864
2.450
5.637
5.981
12.302
RGP CPU (Part I + the rest)
0.219 + 0.350
1.404 + 0.168
5.943 + 2.147
9.975 + 3.049
59.036 + 6.468Table 5
Convergence study on irregular nodes for Example5.2.
Irregular nodes
5
2
175
370
637
976
RGD
r
u
0.0413
0.0062
0.0021
0.0055
0.0175
r
σ
0.1836
0.0399
0.0578
0.0832
0.1138
CPU
0.419
0.831
2.055
2.627
4.387
RGP/RPIM
r
u
0.0386
0.0003
0.0020
0.0038
0.0003
r
σ
0.1758
0.0553
0.0422
0.0610
0.0108
RPIM CPU
1.637
2.486
5.600
6.178
12.814
RGP CPU (Part I + the rest)
0.203
+
0.332
0.832
+
0.734
5.989
+
1.667
10.018
+
2.480
58.798
+
3.672


#### 5.2.4. Discussion on the Computational Results

To further discuss the effectiveness of the RG methods, we discuss the computational results of the cantilever beam problem in this subsection. Figure6 shows a comparison of analytical solution and the present RGD numerical solution for the beam deflection. Two nodal distributions, 91 (13
×
7) regular nodes and 481 (37
×
13) regular nodes, are used. The stress results are also obtained. Figures 8 and 9 illustrate the analytical solution and the RG methods solutions for the stress σ
x and the shear stress τ
x
y of the beam using 481 regular nodes, respectively. Figures 6–8 show the RG methods and RPIM have good performance for the stress σ
x and poor performance for the stress τ
x
y. But in terms of time, the RGD method is better than the RGP method and RPIM method (Tables 4 and 5 and Figure 6). Thus, the RGP method is not considered in the next section.Convergence study with error of displacement: regular nodes (a), irregular nodes (b).
(a)
(b)Displacement13
×
7 (a) and 36
×
13 (b) regular nodes distribution (scaled 400).
(a)
(b)Stressσ
x: analytical results (a), RGP numerical results (b), and RGD numerical results (c).
(a)
(b)
(c)Shear stressτ
x
y: analytical results (a), RGP numerical results (b), and RGD numerical results (c).
(a)
(b)
(c)


#### 5.2.5. A Special Case of the RD Methods

As a special case of the RD methods, the vertices of the Gauss domains coincide with the nodes in regular nodes distribution; see Figure11. Figure 11 shows that a 5
*
5 uniform grid and a 4
×
4 Gauss quadrature in each Gauss domain are used. The Gauss points * in the Gauss domain have the same nodes ∘ (red nodes) with α
=
2. We consider that the process of constructing shape function has three parts: the [
R
(
x
)
P
(
x
)
], solving the inverse of G, and computing [
R
(
x
)
P
(
x
)
]
*
G
-
1. These computational complexities are O
(
(
m
+
n
)
), O
(
(
m
+
n
)
3
), and O
(
(
m
+
n
)
2
), respectively. Table 6 shows the complexity of the integration for a Gauss domain with 16 Gauss points.Table 6
The complexity of the integration for a Gauss domain.
Complexity
Part I
Part II
Part III
Sum
O
(
(
m
+
n
)
)
O
(
(
m
+
n
)
3
)
O
(
(
m
+
n
)
2
)
RPIM
16 times
16 times
16 times
16
*
O
(
(
m
+
n
)
3
)
+
16
*
O
(
(
m
+
n
)
2
)
+
16
*
O
(
m
+
n
)
RG
16 times
1 time
16 times
O
(
(
m
+
n
)
3
)
+
16
*
O
(
(
m
+
n
)
2
)
+
16
*
O
(
m
+
n
)The efficiency ratio is defined asCP
U
RPIM
/
CP
U
RG. To discuss the efficiency ratio, two nodal distributions, 325
(
25
×
13
) regular nodes and 629
(
37
×
17
) regular nodes, are used. Table 7 shows that the efficiency ratios are 11.832 and 13.362, respectively. The efficiency ratios are also expressed as
(23)
efficiency
ratios
=
NG
*
16
*
O
m
+
n
3
+
16
*
O
m
+
n
2
i
i
i
i
i
i
I
i
i
i
i
i
+
16
*
O
m
+
n

m
+
n
3
i
i
+
preprocess
+
postprocess
m
+
n
3
×
NG
*
O
m
+
n
3
+
16
*
O
m
+
n
2
iiiiiiiiiiiiiiiiiiIIiii
+
16
*
O
(
m
+
n
)


m
+
n
3
i
i
i
i
i
i
i
i
i
i
+
preprocess
+
postprocess
m
+
n
3
-
1
.Table 7
Effect of the size of support domain for the RG method (Example4.2).
Grid
25
*
13
37
*
17
r
u
r
σ
CPU
r
u
r
σ
CPU
Schemes
RG
4.7796
e
-
004
0.0204
1.267
7.3559
e
-
005
0.0054
2.373
RPIM
14.991
31.708
Efficiency ratio
11.832
13.362To refine the nodes domain, the NG will be a large number. The elapsed CPU times of compute shape functions is far than the rest. Thus,(24)
efficiency
ratios
⟶
16
*
O
m
+
n
3
i
i
i
i
I
i
i
i
+
16
*
O
m
+
n
2
i
i
i
i
i
I
i
i
i
m
+
n
3
+
16
*
O
m
+
n


m
+
n
3
×
O
m
+
n
3
+
16
*
O
m
+
n
2
i
i
i
i
i
i
I
i
i
m
+
n
3
+
16
*
O
m
+
n
m
+
n
3

-
1
⟶
16
.
The formula of the efficiency ratios fully explains that the ratio is 11.832 for 325
(
25
×
13
) regular nodes (NG
=
288) and 13.362 for 629
(
37
×
17
) regular nodes (NG
=
576). Numerical results show that the rates are 11.832 and 13.362 which gradually close to the theory upper 16 with the increase of NG (the number of the Gauss domains). If we want to improve the accuracy of numerical results with refining the nodes domain, the number of the Gauss points should be added and the efficiency ratios will be improved further.


### 5.3. Plate with a Central Circular Hole

Consider now a plate with a central circular hole subjected to a unidirectional tensile load of 1.0 in thex direction (see Figure 11(a)). This is a typical plane stress problem.The closed-form solution of stresses is(25)
σ
x
=
1
-
a
2
r
2
3
2
cos
⁡
2
θ
+
cos
⁡
4
θ
+
3
a
4
2
r
4
cos
⁡
4
θ
,
σ
y
=
-
a
2
r
2
1
2
cos
⁡
2
θ
-
cos
⁡
4
θ
-
3
a
4
2
r
4
cos
⁡
4
θ
,
τ
x
y
=
-
a
2
r
2
1
2
sin
2
θ
+
sin
4
θ
+
3
a
4
2
r
4
sin
4
θ
,
where (r
,
θ) are polar coordinates, θ is the measured counterclockwise from the positive x-axis, and a is the radius of the hole. The corresponding displacements, in the plain stress case, are given by
(26)
u
r
=
1
4
G
r
κ
-
1
2
+
cos
⁡
2
θ
+
a
2
r
1
+
1
+
κ
cos
⁡
2
θ
iiiiiiiiiiii
-
a
4
r
3
cos
⁡
2
θ
,
u
θ
=
1
4
G
1
-
κ
a
2
r
-
r
-
a
4
r
3
sin
2
θ
,
where G
=
E
/
2
(
1
+
ν
) and κ
=
3
-
ν
/
1
+
ν. The material properties are E
=
1.0
×
1
0
3 and ν
=
0.3.By taking advantage of its symmetry, only a quarter of the model is considered in the analysis (see Figure11(b)). Symmetry conditions are imposed on the left and the bottom edges, and the inner boundary of the hole is traction-free. These boundary conditions include (i) essential boundary conditions on the bottom and left edges on which displacement u
¯ is computed from the exact displacement given in (26) and (ii) natural boundary conditions on the right and top edges on which traction t
¯ is computed from the exact stress given in (25). A typical node distribution with 99 nodes and the corresponding background mesh constructed for the numerical integration are shown in Figure 12(a). For each background cell, a 4
×
4 Gauss quadrature is employed. It is also of great interest to further study the effect of the nodal influence domain size and the computational support domain size. It should be noted that the support domain for all nodes is a circle with varying radius. They are chosen such that the support is small for nodes near the hole and bigger for nodes near the edges. Table 8 presents the relative errors and elapsed CPU times with different parameters α for the 99 nodes. In this case, the computational results show relatively good accuracy and less elapsed CPU time. For convergence studies, four different node distributions with 99 nodes, 289 nodes, 625 nodes, and 1089 nodes are considered in Table 9. The last three node distributions are taken from [23]. Due to singular property for RPIM, α was advised as 1.1 compared with the RGD method. The computational results are presented in Table 9, which also shows that the RGD method possesses high convergence. The elapsed CPU times of the four node distributions fall 9
%, 55
%, 80
%, and 76
%, respectively. The computed displacements and stresses are also shown for 99 nodes’ distribution problem in Figures 10–13. The displacements computed by the RGD method at nodes are plotted and compared to the exact solution (see Figure 12(b)). Figures 13 and 14 compare the analytical solution and RGD solution of stress σ
x and shear stress τ
x
y, respectively. Figure 15(a) compares the stresses σ
x at stations along the left edge of the plate for the RPIM solution, the RGD solution, and the closed-form solution. Figure 15(b) plots the stresses σ
y at stations along the bottom edge of the plate for these methods. From Figures 11–13, it is shown that the RGD has a higher accuracy compared with the RPIM solution. In particular, the RGD solution of stress σ
x at (
0,5
) agrees with the analytical result, but the RPIM result is far away from the analytical result.Table 8
Effect of the size of support domain for the RG methods for the 99 nodes (Example5.3).
α
RGD
RPIM
r
u
r
σ
CPU
r
u
r
σ
CPU
1.1
0.0149
0.1101
0.880
0.0219
0.1005
0.968
1.2
0.0097
0.0881
0.963
0.0205
0.1004
1.052
1.5
0.0094
0.0623
1.175
0.0271
0.1526
1.547
1.8
0.0056
0.0886
1.580
0.0133
0.2052
2.383
2.0
0.0054
0.1154
2.043
0.0115
0.2226
3.052
3.0
0.0044
0.0590
4.600
0.0057
01500
9.286Table 9
Convergence study for Example5.3.
Scheme
99 nodes
289 nodes
625 nodes
1089 nodes
RGD
r
u
0.0149
0.0053
0.0025
0.0016
r
σ
0.1101
0.0534
0.0384
0.0343
CPU
0.880
2.030
3.480
13.516
RPIM
r
u
0.0219
0.0057
0.0031
0.0022
r
σ
0.1005
0.0541
0.0343
0.0283
CPU
0.968
4.523
17.419
55.530Figure 10
The diagram of5
*
5 domains for the special case.Plate with a central circular hole and its model problem: (a) a plate with a hole, (b) a quarter of the plate and boundary conditions.
(a)
(b)Plate with a central circular hole: (a) 99 nodes’ distribution, (b) analytical and RGD solution of displacement.
(a)
(b)Stressσ
x: analytical results (a) and RGD numerical results (b).
(a)
(b)Shear stressτ
x
y: analytical results (a) and RGD numerical results (b).
(a)
(b)Comparisons of RGD and analytical solution: stressσ
x at the stations along the left edge of the plate (a) and stress σ
y at the stations along the bottom edge of the plate (b).
(a)
(b)


## 5.1. Patch Test

The first numerical example is the standard patch test shown in Figure2. The patch test consists of thirteen nodes including five interior nodes. A 2
×
2 rectangular background mesh is used for numerical integration and 4
×
4 Gauss points are used in each mesh. In this patch test, the displacements are prescribed on all outside boundaries by a linear function of x and y on the patches of the dimension. The parameters are taken as E
=
1.0 and ν
=
0.3. The linear displacement functions are u
x
=
0.6
x and u
y
=
0.6
y. RG methods of the patch test require that the displacement of any interior node be given by the same linear function and that the strains and stresses be constant in the patch. The RG methods pass the patch test when linear polynomials are added (m
=
3). However, when polynomial term is not included (m
=
0), the patch test does not easily pass. The same is for the RPIM method; see [21]. Computed results at interior nodes for the RG methods and the RPIM method and the exact results are listed in Table 1.Table 1
Results at interior nodes located irregularly in the 13 nodes’ patch.
Internal node
Coordinates
Exact
RGP/RPIM
RGD
5
(1.0, 1.0)
(
u
x
,
u
y
)
(0.600, 0.600)
(0.601, 0.599)
(0.601, 0.599)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.853, 0.853, 0)
(0.853, 0.853, 0)
10
(0.65, 1.0)
(
u
x
,
u
y
)
(0.390, 0.600)
(0.391, 0.599)
(0.391, 0.599)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.861, 0.853, 0.002)
(0.861, 0.853, 0.002)
11
(0.70, 1.5)
(
u
x
,
u
y
)
(0.420, 0.900)
(0.422, 0.897)
(0.421, 0.897)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.856, 0.855, 0.001)
(0.856, 0.855, 0.001)
12
(1.3, 1.2)
(
u
x
,
u
y
)
(0.780, 0.720)
(0.779, 0.718)
(0.779, 0.718)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.846, 0.852, −0.002)
(0.846, 0.852, −0.002)
13
(1.2, 0.6)
(
u
x
,
u
y
)
(0.720, 0.360)
(0.720, 0.359)
(0.720, 0.359)
(
σ
x
,
σ
y
,
τ
x
y
)
(0.857, 0.857, 0)
(0.854, 0.856, 0)
(0.854, 0.856, 0)Figure 2
A patch with 13 nodes.


## 5.2. Cantilevered Beam

The second example is a cantilever beam problem; see Figure3. Consider a beam of length L and height H subjected to traction at the free end. The beam has a unit thickness, and thus a plane stress problem is considered here. The closed-form solution is available for parabolic traction of force P:
(20)
u
x
=
P
y
6
E
I
6
L
-
3
x
x
+
2
+
ν
y
2
-
H
2
4
,
u
y
=
-
P
6
E
I
3
ν
y
2
L
-
x
+
4
+
5
ν
H
2
x
4
+
3
L
-
x
x
2
,
where the moment of inertia I of the beam is given by I
=
H
3
/
12. The stresses corresponding to the above displacements are
(21)
σ
x
=
P
L
-
x
y
I
,
τ
x
y
=
-
P
2
I
H
2
4
-
y
2
,
σ
y
=
0
.
The beam parameters are taken as E
=
3.0
×
1
0
7,  ν
=
0.3,  H
=
12,  L
=
48, and P
=
1000 in computation.Figure 3
Cantilever beam.To evaluate effects of various parameters, we use the two-norm relative errorsr
u and r
σ for displacement and stress, respectively. They are defined by
(22)
r
u
=
u
-
u
^
2
u
2
,
r
σ
=
σ
^
-
σ
2
σ
2
,
where u
^ and u are the approximate and exact solution of displacements and σ
^ and σ are the approximate and exact value of stresses.


### 5.2.1. Effect of Irregular Node Distribution

Node distributions with 325 irregular nodes and 481 irregular nodes are shown in Figure4. A background mesh is necessary to obtain the numerical integrations of (15). This mesh is independent of nodes for interpolations, while FEM uses the same nodes for both interpolation and numerical integration. For this problem, 10
×
6 and 15
×
6 background cells are used for the 325 irregular nodes’ problem and 481 irregular nodes’ problem, respectively. In each cell, a 4
×
4 Gauss quadrature is used to evaluate the stiffness matrix. Figure 5(a) shows a comparison of the analytical solution and numerical solution of u
y along y
=
0 for this problem with 325 irregular nodes. The plot shows an excellent agreement between the analytical and numerical results for each method when an irregular node distribution is used. Figure 5(b) is the sectional distribution of shear stress τ
x
y along the x
=
24 section for this problem with 325 irregular nodes. The closed-form solution is also plotted for comparison.Meshless models with irregular data points: 325 nodes (a), 481 nodes (b).
(a)
(b)Results for 325 irregular nodes: (a) displacement and (b) shear stress.
(a)
(b)Computational results including error and computational cost are listed in Table2 for the cantilevered beam problem with 325 irregular nodes and 481 irregular nodes. From Figure 5 and Table 2, we can find that the RG methods have better accuracy and are less time-consuming than the RPIM method.Table 2
Computational results for irregular node distribution.
Irregular nodes
RPIM
RGP(Part I + the rest)
RGD
325
r
u
0.0046
0.0013
r
σ
0.0793
0.0582
CPU
3.034
2.917 = 1.424 + 1.493
1.162
481
r
u
0.0080
0.0019
r
σ
0.0824
0.0472
CPU
4.521
6.468 = 4.354 + 2.114
1.668


### 5.2.2. Effect of the Size of Support Domain of RGD

As discussed in Section4.2, the centroid of every Gauss cell has a support domain. The range of the support domain α is an important parameter. The 325 regular nodes’ distribution (25
×
13: 25 nodes in the x direction and 13 nodes in the y direction) is used to study the effect of the size of the support domain for the RG methods.In Table3, we list the numerical results for the RGD method with different α. To compare with the RPIM method, we also list the computational results for the RPIM method with different radiuses of the support domain. From Table 3, we can see that the RGD method has better steadiness, computational efficiency, and computational accuracy than the RPIM method. Particularly, with α
=
4.5, the elapsed CPU time is reduced to 12.7
% (see the data in the boxes) of the previous one. This shows that the size of the support domain has increased and its CPU time takes up a larger proportion of the entire time, in which the consumption surged from 48
% to 76
% of the entire cost. The searching process has become the major expenditure of the entire program. From Table 3, we can find that the RGD method has a stable convergency with the increase of the α, but it does not appear in the RPIM. In this paper, we take the suitable α to be 4.Table 3
Effect of the size of support domain of RGD (Example5.2).
α
RGD
RPIM and RGP
r
u
r
σ
CPU
r
u
r
σ
RPIM CPU
RGP CPU (Part I + the rest)
2.5
0.0513
0.2033
0.413
0.0074
0.0522
1.808
2.151
=
1.023
+
1.128
3.0
0.0129
0.1000
0.646
0.0031
0.0620
3.003
2.537
=
1.131
+
1.406
3.5
0.0109
0.0945
0.733
0.0031
0.0305
4.752
3.277
=
1.798
+
1.479
4.0
0.0006
0.0285
1.182
0.0002
0.0533
7.135
4.071
=
2.608
+
1.463
4.5
0.0004
0.0274
1.336
0.0021
0.0191
10.303
5.195
=
3.720
+
1.475
5.0
0.0002
0.0202
2.112
0.0004
0.0294
14.599
6.451
=
4.911
+
1.540


### 5.2.3. Convergence Study

For convergence studies, five different regular node distributions with 52 (13
×
4: 13 nodes in the x direction and 4 nodes in the y direction) nodes, 175 (25
×
7) nodes, 370 (37
×
10) nodes, 637 (49
×
13) nodes, and 976 (61
×
16) nodes and five different irregular node distributions are considered. The size of the support domain is taken as α
=
4. The computational results are listed in Tables 4 and 5. In Tables 4 and 5, we also list the computational results of the RPIM method. The convergence curves with different node distributions are plotted in Figure 5. In Figure 5, h is the maximum size of node arrangement. From Tables 4 and 5 and Figure 5 we can see that the present method possesses good convergence and less computational times than the RPIM method in general. However, the searching process consumption overtakes 88
% of the entire process in the 976 regular/irregular nodes, which become cumbersome. The oscillatory behavior of the convergence curve of RGD method exists in regular node, and oscillatory behaviors appear in the two node distributions of the RPIM method. Figure 5 shows that RGD method is steadier than the RPIM method. The reason of oscillatory appearance is that the accuracy of most meshless methods is closely related to the integration scheme, the number of Gauss points, the radiuses of the support domain and influence domain, and so on. Thus, it is a difficult task to get the best accuracy for all cases.Table 4
Convergence study on regular nodes for Example5.2.
Regular nodes
13
×
4
25
×
7
37
×
10
49
×
13
61
×
16
RGD
r
u
0.0384
0.0063
0.0035
0.0025
0.0015
r
σ
0.1891
0.0399
0.0676
0.0668
0.0182
CPU
0.440
0.839
1.972
2.730
5.420
RGP/RPIM
r
u
0.0358
0.0028
0.0024
0.0196
0.0240
r
σ
0.2057
0.0522
0.0559
0.1434
0.1738
RPIM CPU
1.864
2.450
5.637
5.981
12.302
RGP CPU (Part I + the rest)
0.219 + 0.350
1.404 + 0.168
5.943 + 2.147
9.975 + 3.049
59.036 + 6.468Table 5
Convergence study on irregular nodes for Example5.2.
Irregular nodes
5
2
175
370
637
976
RGD
r
u
0.0413
0.0062
0.0021
0.0055
0.0175
r
σ
0.1836
0.0399
0.0578
0.0832
0.1138
CPU
0.419
0.831
2.055
2.627
4.387
RGP/RPIM
r
u
0.0386
0.0003
0.0020
0.0038
0.0003
r
σ
0.1758
0.0553
0.0422
0.0610
0.0108
RPIM CPU
1.637
2.486
5.600
6.178
12.814
RGP CPU (Part I + the rest)
0.203
+
0.332
0.832
+
0.734
5.989
+
1.667
10.018
+
2.480
58.798
+
3.672


### 5.2.4. Discussion on the Computational Results

To further discuss the effectiveness of the RG methods, we discuss the computational results of the cantilever beam problem in this subsection. Figure6 shows a comparison of analytical solution and the present RGD numerical solution for the beam deflection. Two nodal distributions, 91 (13
×
7) regular nodes and 481 (37
×
13) regular nodes, are used. The stress results are also obtained. Figures 8 and 9 illustrate the analytical solution and the RG methods solutions for the stress σ
x and the shear stress τ
x
y of the beam using 481 regular nodes, respectively. Figures 6–8 show the RG methods and RPIM have good performance for the stress σ
x and poor performance for the stress τ
x
y. But in terms of time, the RGD method is better than the RGP method and RPIM method (Tables 4 and 5 and Figure 6). Thus, the RGP method is not considered in the next section.Convergence study with error of displacement: regular nodes (a), irregular nodes (b).
(a)
(b)Displacement13
×
7 (a) and 36
×
13 (b) regular nodes distribution (scaled 400).
(a)
(b)Stressσ
x: analytical results (a), RGP numerical results (b), and RGD numerical results (c).
(a)
(b)
(c)Shear stressτ
x
y: analytical results (a), RGP numerical results (b), and RGD numerical results (c).
(a)
(b)
(c)


### 5.2.5. A Special Case of the RD Methods

As a special case of the RD methods, the vertices of the Gauss domains coincide with the nodes in regular nodes distribution; see Figure11. Figure 11 shows that a 5
*
5 uniform grid and a 4
×
4 Gauss quadrature in each Gauss domain are used. The Gauss points * in the Gauss domain have the same nodes ∘ (red nodes) with α
=
2. We consider that the process of constructing shape function has three parts: the [
R
(
x
)
P
(
x
)
], solving the inverse of G, and computing [
R
(
x
)
P
(
x
)
]
*
G
-
1. These computational complexities are O
(
(
m
+
n
)
), O
(
(
m
+
n
)
3
), and O
(
(
m
+
n
)
2
), respectively. Table 6 shows the complexity of the integration for a Gauss domain with 16 Gauss points.Table 6
The complexity of the integration for a Gauss domain.
Complexity
Part I
Part II
Part III
Sum
O
(
(
m
+
n
)
)
O
(
(
m
+
n
)
3
)
O
(
(
m
+
n
)
2
)
RPIM
16 times
16 times
16 times
16
*
O
(
(
m
+
n
)
3
)
+
16
*
O
(
(
m
+
n
)
2
)
+
16
*
O
(
m
+
n
)
RG
16 times
1 time
16 times
O
(
(
m
+
n
)
3
)
+
16
*
O
(
(
m
+
n
)
2
)
+
16
*
O
(
m
+
n
)The efficiency ratio is defined asCP
U
RPIM
/
CP
U
RG. To discuss the efficiency ratio, two nodal distributions, 325
(
25
×
13
) regular nodes and 629
(
37
×
17
) regular nodes, are used. Table 7 shows that the efficiency ratios are 11.832 and 13.362, respectively. The efficiency ratios are also expressed as
(23)
efficiency
ratios
=
NG
*
16
*
O
m
+
n
3
+
16
*
O
m
+
n
2
i
i
i
i
i
i
I
i
i
i
i
i
+
16
*
O
m
+
n

m
+
n
3
i
i
+
preprocess
+
postprocess
m
+
n
3
×
NG
*
O
m
+
n
3
+
16
*
O
m
+
n
2
iiiiiiiiiiiiiiiiiiIIiii
+
16
*
O
(
m
+
n
)


m
+
n
3
i
i
i
i
i
i
i
i
i
i
+
preprocess
+
postprocess
m
+
n
3
-
1
.Table 7
Effect of the size of support domain for the RG method (Example4.2).
Grid
25
*
13
37
*
17
r
u
r
σ
CPU
r
u
r
σ
CPU
Schemes
RG
4.7796
e
-
004
0.0204
1.267
7.3559
e
-
005
0.0054
2.373
RPIM
14.991
31.708
Efficiency ratio
11.832
13.362To refine the nodes domain, the NG will be a large number. The elapsed CPU times of compute shape functions is far than the rest. Thus,(24)
efficiency
ratios
⟶
16
*
O
m
+
n
3
i
i
i
i
I
i
i
i
+
16
*
O
m
+
n
2
i
i
i
i
i
I
i
i
i
m
+
n
3
+
16
*
O
m
+
n


m
+
n
3
×
O
m
+
n
3
+
16
*
O
m
+
n
2
i
i
i
i
i
i
I
i
i
m
+
n
3
+
16
*
O
m
+
n
m
+
n
3

-
1
⟶
16
.
The formula of the efficiency ratios fully explains that the ratio is 11.832 for 325
(
25
×
13
) regular nodes (NG
=
288) and 13.362 for 629
(
37
×
17
) regular nodes (NG
=
576). Numerical results show that the rates are 11.832 and 13.362 which gradually close to the theory upper 16 with the increase of NG (the number of the Gauss domains). If we want to improve the accuracy of numerical results with refining the nodes domain, the number of the Gauss points should be added and the efficiency ratios will be improved further.


## 5.2.1. Effect of Irregular Node Distribution

Node distributions with 325 irregular nodes and 481 irregular nodes are shown in Figure4. A background mesh is necessary to obtain the numerical integrations of (15). This mesh is independent of nodes for interpolations, while FEM uses the same nodes for both interpolation and numerical integration. For this problem, 10
×
6 and 15
×
6 background cells are used for the 325 irregular nodes’ problem and 481 irregular nodes’ problem, respectively. In each cell, a 4
×
4 Gauss quadrature is used to evaluate the stiffness matrix. Figure 5(a) shows a comparison of the analytical solution and numerical solution of u
y along y
=
0 for this problem with 325 irregular nodes. The plot shows an excellent agreement between the analytical and numerical results for each method when an irregular node distribution is used. Figure 5(b) is the sectional distribution of shear stress τ
x
y along the x
=
24 section for this problem with 325 irregular nodes. The closed-form solution is also plotted for comparison.Meshless models with irregular data points: 325 nodes (a), 481 nodes (b).
(a)
(b)Results for 325 irregular nodes: (a) displacement and (b) shear stress.
(a)
(b)Computational results including error and computational cost are listed in Table2 for the cantilevered beam problem with 325 irregular nodes and 481 irregular nodes. From Figure 5 and Table 2, we can find that the RG methods have better accuracy and are less time-consuming than the RPIM method.Table 2
Computational results for irregular node distribution.
Irregular nodes
RPIM
RGP(Part I + the rest)
RGD
325
r
u
0.0046
0.0013
r
σ
0.0793
0.0582
CPU
3.034
2.917 = 1.424 + 1.493
1.162
481
r
u
0.0080
0.0019
r
σ
0.0824
0.0472
CPU
4.521
6.468 = 4.354 + 2.114
1.668


## 5.2.2. Effect of the Size of Support Domain of RGD

As discussed in Section4.2, the centroid of every Gauss cell has a support domain. The range of the support domain α is an important parameter. The 325 regular nodes’ distribution (25
×
13: 25 nodes in the x direction and 13 nodes in the y direction) is used to study the effect of the size of the support domain for the RG methods.In Table3, we list the numerical results for the RGD method with different α. To compare with the RPIM method, we also list the computational results for the RPIM method with different radiuses of the support domain. From Table 3, we can see that the RGD method has better steadiness, computational efficiency, and computational accuracy than the RPIM method. Particularly, with α
=
4.5, the elapsed CPU time is reduced to 12.7
% (see the data in the boxes) of the previous one. This shows that the size of the support domain has increased and its CPU time takes up a larger proportion of the entire time, in which the consumption surged from 48
% to 76
% of the entire cost. The searching process has become the major expenditure of the entire program. From Table 3, we can find that the RGD method has a stable convergency with the increase of the α, but it does not appear in the RPIM. In this paper, we take the suitable α to be 4.Table 3
Effect of the size of support domain of RGD (Example5.2).
α
RGD
RPIM and RGP
r
u
r
σ
CPU
r
u
r
σ
RPIM CPU
RGP CPU (Part I + the rest)
2.5
0.0513
0.2033
0.413
0.0074
0.0522
1.808
2.151
=
1.023
+
1.128
3.0
0.0129
0.1000
0.646
0.0031
0.0620
3.003
2.537
=
1.131
+
1.406
3.5
0.0109
0.0945
0.733
0.0031
0.0305
4.752
3.277
=
1.798
+
1.479
4.0
0.0006
0.0285
1.182
0.0002
0.0533
7.135
4.071
=
2.608
+
1.463
4.5
0.0004
0.0274
1.336
0.0021
0.0191
10.303
5.195
=
3.720
+
1.475
5.0
0.0002
0.0202
2.112
0.0004
0.0294
14.599
6.451
=
4.911
+
1.540


## 5.2.3. Convergence Study

For convergence studies, five different regular node distributions with 52 (13
×
4: 13 nodes in the x direction and 4 nodes in the y direction) nodes, 175 (25
×
7) nodes, 370 (37
×
10) nodes, 637 (49
×
13) nodes, and 976 (61
×
16) nodes and five different irregular node distributions are considered. The size of the support domain is taken as α
=
4. The computational results are listed in Tables 4 and 5. In Tables 4 and 5, we also list the computational results of the RPIM method. The convergence curves with different node distributions are plotted in Figure 5. In Figure 5, h is the maximum size of node arrangement. From Tables 4 and 5 and Figure 5 we can see that the present method possesses good convergence and less computational times than the RPIM method in general. However, the searching process consumption overtakes 88
% of the entire process in the 976 regular/irregular nodes, which become cumbersome. The oscillatory behavior of the convergence curve of RGD method exists in regular node, and oscillatory behaviors appear in the two node distributions of the RPIM method. Figure 5 shows that RGD method is steadier than the RPIM method. The reason of oscillatory appearance is that the accuracy of most meshless methods is closely related to the integration scheme, the number of Gauss points, the radiuses of the support domain and influence domain, and so on. Thus, it is a difficult task to get the best accuracy for all cases.Table 4
Convergence study on regular nodes for Example5.2.
Regular nodes
13
×
4
25
×
7
37
×
10
49
×
13
61
×
16
RGD
r
u
0.0384
0.0063
0.0035
0.0025
0.0015
r
σ
0.1891
0.0399
0.0676
0.0668
0.0182
CPU
0.440
0.839
1.972
2.730
5.420
RGP/RPIM
r
u
0.0358
0.0028
0.0024
0.0196
0.0240
r
σ
0.2057
0.0522
0.0559
0.1434
0.1738
RPIM CPU
1.864
2.450
5.637
5.981
12.302
RGP CPU (Part I + the rest)
0.219 + 0.350
1.404 + 0.168
5.943 + 2.147
9.975 + 3.049
59.036 + 6.468Table 5
Convergence study on irregular nodes for Example5.2.
Irregular nodes
5
2
175
370
637
976
RGD
r
u
0.0413
0.0062
0.0021
0.0055
0.0175
r
σ
0.1836
0.0399
0.0578
0.0832
0.1138
CPU
0.419
0.831
2.055
2.627
4.387
RGP/RPIM
r
u
0.0386
0.0003
0.0020
0.0038
0.0003
r
σ
0.1758
0.0553
0.0422
0.0610
0.0108
RPIM CPU
1.637
2.486
5.600
6.178
12.814
RGP CPU (Part I + the rest)
0.203
+
0.332
0.832
+
0.734
5.989
+
1.667
10.018
+
2.480
58.798
+
3.672


## 5.2.4. Discussion on the Computational Results

To further discuss the effectiveness of the RG methods, we discuss the computational results of the cantilever beam problem in this subsection. Figure6 shows a comparison of analytical solution and the present RGD numerical solution for the beam deflection. Two nodal distributions, 91 (13
×
7) regular nodes and 481 (37
×
13) regular nodes, are used. The stress results are also obtained. Figures 8 and 9 illustrate the analytical solution and the RG methods solutions for the stress σ
x and the shear stress τ
x
y of the beam using 481 regular nodes, respectively. Figures 6–8 show the RG methods and RPIM have good performance for the stress σ
x and poor performance for the stress τ
x
y. But in terms of time, the RGD method is better than the RGP method and RPIM method (Tables 4 and 5 and Figure 6). Thus, the RGP method is not considered in the next section.Convergence study with error of displacement: regular nodes (a), irregular nodes (b).
(a)
(b)Displacement13
×
7 (a) and 36
×
13 (b) regular nodes distribution (scaled 400).
(a)
(b)Stressσ
x: analytical results (a), RGP numerical results (b), and RGD numerical results (c).
(a)
(b)
(c)Shear stressτ
x
y: analytical results (a), RGP numerical results (b), and RGD numerical results (c).
(a)
(b)
(c)


## 5.2.5. A Special Case of the RD Methods

As a special case of the RD methods, the vertices of the Gauss domains coincide with the nodes in regular nodes distribution; see Figure11. Figure 11 shows that a 5
*
5 uniform grid and a 4
×
4 Gauss quadrature in each Gauss domain are used. The Gauss points * in the Gauss domain have the same nodes ∘ (red nodes) with α
=
2. We consider that the process of constructing shape function has three parts: the [
R
(
x
)
P
(
x
)
], solving the inverse of G, and computing [
R
(
x
)
P
(
x
)
]
*
G
-
1. These computational complexities are O
(
(
m
+
n
)
), O
(
(
m
+
n
)
3
), and O
(
(
m
+
n
)
2
), respectively. Table 6 shows the complexity of the integration for a Gauss domain with 16 Gauss points.Table 6
The complexity of the integration for a Gauss domain.
Complexity
Part I
Part II
Part III
Sum
O
(
(
m
+
n
)
)
O
(
(
m
+
n
)
3
)
O
(
(
m
+
n
)
2
)
RPIM
16 times
16 times
16 times
16
*
O
(
(
m
+
n
)
3
)
+
16
*
O
(
(
m
+
n
)
2
)
+
16
*
O
(
m
+
n
)
RG
16 times
1 time
16 times
O
(
(
m
+
n
)
3
)
+
16
*
O
(
(
m
+
n
)
2
)
+
16
*
O
(
m
+
n
)The efficiency ratio is defined asCP
U
RPIM
/
CP
U
RG. To discuss the efficiency ratio, two nodal distributions, 325
(
25
×
13
) regular nodes and 629
(
37
×
17
) regular nodes, are used. Table 7 shows that the efficiency ratios are 11.832 and 13.362, respectively. The efficiency ratios are also expressed as
(23)
efficiency
ratios
=
NG
*
16
*
O
m
+
n
3
+
16
*
O
m
+
n
2
i
i
i
i
i
i
I
i
i
i
i
i
+
16
*
O
m
+
n

m
+
n
3
i
i
+
preprocess
+
postprocess
m
+
n
3
×
NG
*
O
m
+
n
3
+
16
*
O
m
+
n
2
iiiiiiiiiiiiiiiiiiIIiii
+
16
*
O
(
m
+
n
)


m
+
n
3
i
i
i
i
i
i
i
i
i
i
+
preprocess
+
postprocess
m
+
n
3
-
1
.Table 7
Effect of the size of support domain for the RG method (Example4.2).
Grid
25
*
13
37
*
17
r
u
r
σ
CPU
r
u
r
σ
CPU
Schemes
RG
4.7796
e
-
004
0.0204
1.267
7.3559
e
-
005
0.0054
2.373
RPIM
14.991
31.708
Efficiency ratio
11.832
13.362To refine the nodes domain, the NG will be a large number. The elapsed CPU times of compute shape functions is far than the rest. Thus,(24)
efficiency
ratios
⟶
16
*
O
m
+
n
3
i
i
i
i
I
i
i
i
+
16
*
O
m
+
n
2
i
i
i
i
i
I
i
i
i
m
+
n
3
+
16
*
O
m
+
n


m
+
n
3
×
O
m
+
n
3
+
16
*
O
m
+
n
2
i
i
i
i
i
i
I
i
i
m
+
n
3
+
16
*
O
m
+
n
m
+
n
3

-
1
⟶
16
.
The formula of the efficiency ratios fully explains that the ratio is 11.832 for 325
(
25
×
13
) regular nodes (NG
=
288) and 13.362 for 629
(
37
×
17
) regular nodes (NG
=
576). Numerical results show that the rates are 11.832 and 13.362 which gradually close to the theory upper 16 with the increase of NG (the number of the Gauss domains). If we want to improve the accuracy of numerical results with refining the nodes domain, the number of the Gauss points should be added and the efficiency ratios will be improved further.


## 5.3. Plate with a Central Circular Hole

Consider now a plate with a central circular hole subjected to a unidirectional tensile load of 1.0 in thex direction (see Figure 11(a)). This is a typical plane stress problem.The closed-form solution of stresses is(25)
σ
x
=
1
-
a
2
r
2
3
2
cos
⁡
2
θ
+
cos
⁡
4
θ
+
3
a
4
2
r
4
cos
⁡
4
θ
,
σ
y
=
-
a
2
r
2
1
2
cos
⁡
2
θ
-
cos
⁡
4
θ
-
3
a
4
2
r
4
cos
⁡
4
θ
,
τ
x
y
=
-
a
2
r
2
1
2
sin
2
θ
+
sin
4
θ
+
3
a
4
2
r
4
sin
4
θ
,
where (r
,
θ) are polar coordinates, θ is the measured counterclockwise from the positive x-axis, and a is the radius of the hole. The corresponding displacements, in the plain stress case, are given by
(26)
u
r
=
1
4
G
r
κ
-
1
2
+
cos
⁡
2
θ
+
a
2
r
1
+
1
+
κ
cos
⁡
2
θ
iiiiiiiiiiii
-
a
4
r
3
cos
⁡
2
θ
,
u
θ
=
1
4
G
1
-
κ
a
2
r
-
r
-
a
4
r
3
sin
2
θ
,
where G
=
E
/
2
(
1
+
ν
) and κ
=
3
-
ν
/
1
+
ν. The material properties are E
=
1.0
×
1
0
3 and ν
=
0.3.By taking advantage of its symmetry, only a quarter of the model is considered in the analysis (see Figure11(b)). Symmetry conditions are imposed on the left and the bottom edges, and the inner boundary of the hole is traction-free. These boundary conditions include (i) essential boundary conditions on the bottom and left edges on which displacement u
¯ is computed from the exact displacement given in (26) and (ii) natural boundary conditions on the right and top edges on which traction t
¯ is computed from the exact stress given in (25). A typical node distribution with 99 nodes and the corresponding background mesh constructed for the numerical integration are shown in Figure 12(a). For each background cell, a 4
×
4 Gauss quadrature is employed. It is also of great interest to further study the effect of the nodal influence domain size and the computational support domain size. It should be noted that the support domain for all nodes is a circle with varying radius. They are chosen such that the support is small for nodes near the hole and bigger for nodes near the edges. Table 8 presents the relative errors and elapsed CPU times with different parameters α for the 99 nodes. In this case, the computational results show relatively good accuracy and less elapsed CPU time. For convergence studies, four different node distributions with 99 nodes, 289 nodes, 625 nodes, and 1089 nodes are considered in Table 9. The last three node distributions are taken from [23]. Due to singular property for RPIM, α was advised as 1.1 compared with the RGD method. The computational results are presented in Table 9, which also shows that the RGD method possesses high convergence. The elapsed CPU times of the four node distributions fall 9
%, 55
%, 80
%, and 76
%, respectively. The computed displacements and stresses are also shown for 99 nodes’ distribution problem in Figures 10–13. The displacements computed by the RGD method at nodes are plotted and compared to the exact solution (see Figure 12(b)). Figures 13 and 14 compare the analytical solution and RGD solution of stress σ
x and shear stress τ
x
y, respectively. Figure 15(a) compares the stresses σ
x at stations along the left edge of the plate for the RPIM solution, the RGD solution, and the closed-form solution. Figure 15(b) plots the stresses σ
y at stations along the bottom edge of the plate for these methods. From Figures 11–13, it is shown that the RGD has a higher accuracy compared with the RPIM solution. In particular, the RGD solution of stress σ
x at (
0,5
) agrees with the analytical result, but the RPIM result is far away from the analytical result.Table 8
Effect of the size of support domain for the RG methods for the 99 nodes (Example5.3).
α
RGD
RPIM
r
u
r
σ
CPU
r
u
r
σ
CPU
1.1
0.0149
0.1101
0.880
0.0219
0.1005
0.968
1.2
0.0097
0.0881
0.963
0.0205
0.1004
1.052
1.5
0.0094
0.0623
1.175
0.0271
0.1526
1.547
1.8
0.0056
0.0886
1.580
0.0133
0.2052
2.383
2.0
0.0054
0.1154
2.043
0.0115
0.2226
3.052
3.0
0.0044
0.0590
4.600
0.0057
01500
9.286Table 9
Convergence study for Example5.3.
Scheme
99 nodes
289 nodes
625 nodes
1089 nodes
RGD
r
u
0.0149
0.0053
0.0025
0.0016
r
σ
0.1101
0.0534
0.0384
0.0343
CPU
0.880
2.030
3.480
13.516
RPIM
r
u
0.0219
0.0057
0.0031
0.0022
r
σ
0.1005
0.0541
0.0343
0.0283
CPU
0.968
4.523
17.419
55.530Figure 10
The diagram of5
*
5 domains for the special case.Plate with a central circular hole and its model problem: (a) a plate with a hole, (b) a quarter of the plate and boundary conditions.
(a)
(b)Plate with a central circular hole: (a) 99 nodes’ distribution, (b) analytical and RGD solution of displacement.
(a)
(b)Stressσ
x: analytical results (a) and RGD numerical results (b).
(a)
(b)Shear stressτ
x
y: analytical results (a) and RGD numerical results (b).
(a)
(b)Comparisons of RGD and analytical solution: stressσ
x at the stations along the left edge of the plate (a) and stress σ
y at the stations along the bottom edge of the plate (b).
(a)
(b)


## 6. Conclusions

In this paper, we have studied researching Gauss points and reordering Gauss domain technique (GD methods), where the computational processes are more reasonable. Secondly, this technique is compatible with most of the common RPIM methods, and this aspect is worth studying. Finally, 3D problems could be solved with this scheme to save the computation time.



---

*Source: 219538-2014-12-23.xml*