# Error Estimate and Adaptive Refinement in Mixed Discrete Least Squares Meshless Method

Authors: J. Amani; A. Saboor Bagherzadeh; T. Rabczuk
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The node moving and multistage node enrichment adaptive refinement procedures are extended in mixed discrete least squares meshless (MDLSM) method for efficient analysis of elasticity problems. In the formulation of MDLSM method, mixed formulation is accepted to avoid second-order differentiation of shape functions and to obtain displacements and stresses simultaneously. In the refinement procedures, a robust error estimator based on the value of the least square residuals functional of the governing differential equations and its boundaries at nodal points is used which is inherently available from the MDLSM formulation and can efficiently identify the zones with higher numerical errors. The results are compared with the refinement procedures in the irreducible formulation of discrete least squares meshless (DLSM) method and show the accuracy and efficiency of the proposed procedures. Also, the comparison of the error norms and convergence rate show the fidelity of the proposed adaptive refinement procedures in the MDLSM method.


---

## Body

## 1. Introduction

Adaptivity needs appropriate numerical solutions of problems in order to describe high gradient regions and an anisotropic behavior of the solution. In an adaptive procedure, a good error estimator plays a very important role. The error estimation in numerical methods is obviously as old as the numerical computations themselves. The earliest paper by Richardson [1] proposed an error estimation procedure to use in finite difference method. This was followed in finite element method (FEM) [2]. A particular strength of the FEM is the well-developed theories of error estimation and adaptivity. Three h-refinement procedures, namely, mesh movement, mesh enrichment, and remeshing have been proposed for adaptivity [3]. In the mesh movement, the total number of nodes remains constant, but the location of the nodes can change in order to achieve a better overall distribution of the error. In the mesh enrichment, the original nodes hold fix and hierarchical nodes or simply more nodes add to the problem domain based on error distribution. In the remeshing, a completely new nodes is constructed based on the information acquired from the previous computation, and hence, it is required to implement a suitable node generator. on one hand, the mesh movement is more suitable than mesh enrichment because the problem scale remains constant, and on the other hand, its interpolations become too distort in the mesh-based methods [4].In order to avoid these problems, an alternative approach, known as meshless methods (MMs), has been developed in recent decades to discretize a continuum body only by a finite number of nodes. In MMs the unknowns are interpolated from the nodal values that constitute the problem degrees-of-freedom. The main advantage of MMs is the fact that the interpolation accuracy is much less affected by the nodal distribution. Many meshless methods have been introduced since Gingold and Monaghan [5] proposed smoothed particle hydrodynamics (SPH) method. Nayroles et al. [6] implemented the diffuse element method (DEM). Belytschko et al. [7] presented the Element-Free Galerkin (EFG) method. Liu et al. [8] suggested the reproducing kernel particle method (RKPM). The other meshless methods that have been developed in recent years are the Finite Point (FP) method [9], the HP clouds method [10], the meshless local Petrov-Galerkin (MLPG) method [11], the local boundary integral equation (LBIE) method [12], the finite cloud (FC) method [13] and the discrete least squares meshless (DLSM) method [14].Researchers used the advantages of MMs for developing efficient error estimate and adaptivity procedures. Rabczuk and Belytschko [15, 16] proposed an adaptive continuum/discrete crack approach for meshfree particle methods and also an adaptivity procedure for structured meshfree particle methods in 2D and 3D problems. Yoon et al. [17] worked on enriched meshfree collocation method with diffuse derivatives for elastic fracture. Zi et al. [18] investigated extended meshfree methods without branch enrichment for cohesive cracks. Bordas et al. [19] proposed enriched meshfree methods without asymptotic enrichment for 3D nonlinear fracture mechanics. Rabczuk and Samaniego [20] worked on discontinuous modelling of shear bands using adaptive meshfree methods. Zhuang et al. [21–23] investigated error control in the EFG method and adaptivity for structured meshfree particle methods in 2D and 3D problems. The DLSM method was extendted for error estimate and adaptivity in solid [24, 25] and fluid [26] problems.Two different formulations, namely irreducible and mixed formulations have been introduced and used for the solution of engineering problems. With the mixed formulation, the continuity requirement decreases by one order compared to the irreducible formulation [27]. Use of mixed formula may result in an improved approximation, in particular, for the gradient variables, which in turn could result in higher accuracy than possible with the irreducible formulation [27]. In the standard mixed FEM, in order to obtain a coefficient matrix which leads to the system of equations with a unique and stable solution, the polynomial functions chosen for approximation of stresses and displacements must satisfy the Ladyzhenskaya-Babuška-Brezzi (LBB or inf-sup) condition [28, 29]. The stability of mixed discretization does not allow FEM to choose independently the approximation spaces, so these spaces are restricted in the stability condition which is known as the LBB condition. However, the least squares approximation has the advantage that it does not require satisfying the LBB condition [30–32]. Hence, this advantage was used by Amani et al. [33] to implement a mixed meshless method named mixed discrete least squares Meshless (MDLSM) method which is formulated based on the least squares residuals functional of the governing partial differential equations of planar elasticity problem and its boundary conditions at the nodal points, and hence, it is stable and is not required to satisfy the LBB condition between the displacements and stresses approximations. Hence, the approximation spaces of the displacements and stresses can be choosen independently while they are obtained simuletanously.In this paper, the MDLSM method is extended for the residual based error estimation and for the two types of adaptive refinement procedures. The node moving adaptive refinement procedure based on the spring analogy [24] and the node enrichment adaptive refinement procedure [25] are formulated and used in the MDLSM method for efficient analysis of the elasticity problems.The present paper is organized as follows. Formulation of the mixed discrete least squares meshless method for solving the planar elasticity problems is given in Section2. In Section 3, an error estimator based on the least square functional residuals is formulated for the MDLSM method to use in the node moving and node enrichment adaptive refinement procedures. In Section 4, we present some numerical benchmark examples which illustrate the proposed adaptive refinement process as well as the efficiency of the error estimator. Finally, some concluding remarks are addressed in Section 5.


## 2. Formulation of Mixed Discrete Least Squares Meshless Method for Elasticity

Consider the following two-dimensional linear elasticity problem(1)
-
μ
Δ
u
+
(
λ
+
μ
)
∇
(
∇
·
u
)
=
f
in
Ω
,
with displacement and traction boundary conditions as follow:
(2)
u
=
u
-
,
v
=
v
-
,
in
Γ
u
,
σ
x
n
x
+
τ
x
y
n
y
=
t
-
x
,
τ
x
y
n
x
+
σ
y
n
y
=
t
-
y
,
aaaaaaaaaaaaaaaaaaaaaaaaaa
in
Γ
t
,
where Ω is a bounded domain representing the region occupied by an elastic body, and λ, μ are the Lame constants which are defined as
(3)
μ
=
E
2
(
1
+
ν
)
>
0
,
λ
=
E
ν
(
1
-
2
ν
)
(
1
+
ν
)
>
0
,
where ν is the Poisson ratio, E is the Young modulus, and Γ
u, Γ
t are the displacement and traction boundaries, respectively. u
-, v
-, t
-
x, and t
-
y prescribed respectively the displacements and tractions in the x and y directions and n
x, n
y are direction cosines of the normal vector to the boundary.By using the following definition of stresses in terms of the displacement components:(4)
σ
x
=
(
λ
+
2
μ
)
∂
u
∂
x
+
λ
∂
v
∂
y
,
σ
y
=
λ
∂
u
∂
x
+
(
λ
+
2
μ
)
∂
v
∂
y
,
τ
x
y
=
μ
(
∂
u
∂
y
+
∂
v
∂
x
)
,
we can rewrite (1) in term of stresses as
(5)
∂
σ
x
∂
x
+
∂
τ
x
y
∂
y
=
-
f
x
in
Ω
,
∂
τ
x
y
∂
x
+
∂
σ
y
∂
y
=
-
f
y
in
Ω
.The compact form of (1) can be written by substituting (4)-(5) into the second-order problem of (1) in the form of
(6)
L
(
ϕ
)
+
f
=
0
,
where L
(
·
) is a first-order differential operator defined as
(7)
L
(
·
)
=
L
1
(
·
)
x
+
L
2
(
·
)
y
+
L
3
(
·
)
,
and ϕ is the vector of unknowns defined as
(8)
ϕ
=
[
u
v
σ
x
σ
y
τ
x
y
]
T
,
and vector f contains the forcing terms which has the form
(9)
f
=
[
0
0
0
-
f
x
-
f
y
]
T
.
In (7), L
1, L
2, and L
3 are defined by the following matrices:
(10)
L
1
=
(
λ
+
2
μ
0
0
0
0
λ
0
0
0
0
0
μ
0
0
0
0
0
1
0
0
0
0
0
0
1
)
,
L
2
=
(
0
λ
0
0
0
0
λ
+
2
μ
0
0
0
μ
0
0
0
0
0
0
0
0
1
0
0
0
1
0
)
,
L
3
=
(
0
0
-
1
0
0
0
0
0
-
1
0
0
0
0
0
-
1
0
0
0
0
0
0
0
0
0
0
)
.
The displacement and traction boundary conditions (2) can be written in terms of the unknown vector ϕ as
(11)
D
ϕ
-
f
-
=
0
,
where D and f
- are defined as follow:
(12)
D
=
(
1
0
0
0
0
0
1
0
0
0
0
0
n
x
0
n
y
0
0
0
n
y
n
x
)
,
f
-
=
(
u
-
v
-
t
-
x
t
-
y
)
.
The plane elasticity problem is now defined as solving the first-order differential equation subjected only to the Dirichlet type boundary condition
(13)
L
1
(
ϕ
)
x
+
L
2
(
ϕ
)
y
+
L
3
(
ϕ
)
+
f
=
0
,
in
Ω
,
D
ϕ
-
f
-
=
0
,
on
⁡
Γ
.The application of the proposed MDLSM method for solving problem of (13) starts with the definition of residuals as follows:
(14)
ℛ
Ω
=
L
(
ϕ
)
+
f
in
Ω
,
ℛ
Γ
=
D
ϕ
-
f
-
on
⁡
Γ
,
where ℛ
Ω and ℛ
Γ are domain and boundary residuals, respectively.Now the penalty approach is used to form the least square residuals functional which is defined as(15)
I
=
∑
k
=
1
M
d
ℛ
Ω
T
ℛ
Ω
+
α
∑
k
=
1
M
b
ℛ
Γ
T
ℛ
Γ
,
M
=
M
d
+
M
b
,
where M is the total number of sampling (or collocation) points, M
d is total number of domain sampling points, M
b is total number of boundary sampling points, and the penalty coefficient α is a positive scalar constant that must be large enough in order to impose the essential boundary condition with the desired accuracy. A note should be made here regarding the value of the penalty parameters. To impose the boundary conditions exactly, the penalty factor must be infinite, which is not possible in practical numerical analysis. Therefore, the boundary conditions could not be satisfied exactly but only approximately. In general, the use of a larger penalty factor will lead to better enforcement of the constraint. The proper value of the penalty parameter is determined prior to the main calculation via a trial and error process and it is problem dependent.Minimizing the functional in (15) with respect to the nodal unknown vector ϕ leads to the following system of equation:
(16)
K
Φ
=
F
,
where
(17)
K
i
j
=
∑
ℓ
=
1
M
d
‍
[
L
(
N
i
)
]
ℓ
T
[
L
(
N
j
)
]
ℓ
+
α
∑
ℓ
=
1
M
b
‍
[
D
(
N
i
)
]
ℓ
T
[
D
(
N
j
)
]
ℓ
,
F
i
=
∑
ℓ
=
1
M
d
‍
[
L
(
N
i
)
]
ℓ
T
f
ℓ
+
α
∑
ℓ
=
1
M
b
‍
[
D
(
N
i
)
]
ℓ
T
f
-
ℓ
,
and Φ is unknown matrix that contains displacements and stresses of all nodes. F is the right hand side vector and the stiffness matrix K in (16) is n
DOF
×
n
DOF square matrix where n
DOF is the number of unknowns per each node and N is the moving least squares (MLS) shape functions. The proposed MDLSM method has 2.5 times unknowns compared to the irreducible DLSM method with only displacements as unknowns. This should drastically reduce the computational efficiency by (
2.5
)
β (where β depends on the type of linear solver used) times. But, since matrix K is symmetric and positive definite, therefore, the final system of equation can be solved directly via efficient solvers. The MDLSM formulation for the plane elasticity problem has the following advantages that increase its efficiency [33].(1)
The order of shape function derivatives is reduced by one order, thus complex and costly second-order derivative calculations of the MLS shape function in the irreducible DLSM method are avoided.
(2)
The stresses are obtained directly, while calculation of stresses in the irreducible DLSM method requires some postprocessing.
(3)
Only a linear complete polynomial basis is needed to construct the MLS shape functions, while in the irreducible DLSM formulation, the second order polynomial basis is required to achieve second-order consistency due to the presence of second-order derivatives in the irreducible formulation. This leads to lower computational effort for mixed method in construction of MLS shape functions and its derivatives compared to irreducible one. Furthermore, this in turn leads to the smaller number of points to be included in the support domain for the construction of the MLS shape functions in the mixed method and hence, less computational effort.
(4)
Both the displacement and stress boundary conditions are of the Dirichlet type which requires the specification of a single penalty parameter if a penalty method is used to enforce them as used in this work. In the irreducible DLSM method both Dirichlet and Neumann type boundary conditions are required for the determination of two types of penalty coefficient.
(5)
The mixed formulation, when used with the standard weighted residual methods both mesh-based and meshless forms, requires the LBB condition because the resulting problem is a saddle point problem. The least squares method, however, is a minimization method and therefore is not subject to the LBB condition.


## 3. Error Estimator and Adaptive Refinement

### 3.1. Error Estimator

In numerical methods, a problem is solved by discretization of the problem domain into the subdomains, hence, the governing equations only apply into these subdomains, so numerical methods always come with discretization error. Discretization error is one of the most important challenges in the numerical methods. Discretization error is theoretically decreased by refining the discretization domain but perfunctory refinement imposes the heavy computational cost without supplying the expected accuracy. Adaptive refinement methods mean balances between refinement procedure and its computational cost. These methods only refine locally the regions of the domain which has higher error. Adaptive procedure has two main parts: error estimation and adaptive refinement. Any success adaptive refinement needs a reliable error estimation procedure. Real error distribution can not be practicably used because the exact solution is not available for any practical problems. Several methods are used for error estimation with different numerical methods and these methods are categorized into two classes, namely, the residual based method [34] and recovery based method [35]. In residual based method, the residuals of differential equation and its boundaries are used as a criterion of error. The gradient of the solutions is used in recovery based method as the error criterion.In this paper, the relative least square residuals functional for each node is defined as follows:(18)
e
=
I
(
U
t
T
U
t
)
,
where I is the least square residuals functional in (15) and U
t is unknowns obtained from the main solution. It is noticed that most of the computations of the least square residuals functional can be obtained from the main solution of the MDLSM method.


### 3.2. Adaptive Refinement Procedures

#### 3.2.1. Node Moving

Mesh movement strategy can be easily and efficiently used with meshless methods since no element distortion is associated with the method. It should be noted that the mesh movement technique can be used in conjunction with the MDLSM method to adaptively adjust nodal points to improve the quality of the solution obtained with a prespecified number of nodal points. Here, a nodal refinement procedure is used that is called node moving adaptive refinement approach. When a node refinement is required, springs of prescribed stiffness are placed between each pair of nodes belonging to the same subdomain and the nodes are then moved until the spring system is in equilibrium.In the node moving procedure, first, all nodes are connected with springs in which the neighbor nodes are defined using Voronoi diagram [36] (see Figure 1). Voronoi diagram is defined as
(19)
T
i
=
{
x
∈
ℝ
:
d
(
x
,
x
i
)
<
d
(
x
,
x
j
)
,
j
≠
i
}
,
where d
(
x
,
x
i
) is Euclidean distance between x and x
i. The above equation means that neighbor nodes to node i are the closest nodes to the node i rather than other nodes.Figure 1
Voronoi diagram of neighboring nodes and spring connections.Spring forces are defined as(20)
b
i
j
=
c
i
j
(
x
i
-
x
j
)
,
where c
i
j is stiffness of spring between i, j and x
i, x
j are coordinates of i, j in equilibrium, respectively. The free body diagram is shown in Figure 2.Figure 2
Free body diagram between two nodesi and j.Spring stiffness is defined as a function of errors between two pointsi, j as follows:
(21)
c
i
j
=
(
e
i
+
e
j
)
d
i
j
,
where e
i and e
j are the values of the error estimators obtained from (18) at nodes i and j, respectively, and d
i
j is distance between these two connected nodes. In matrix form we have
(22)
(
c
i
j
0
-
c
i
j
0
0
c
i
j
0
-
c
i
j
-
c
i
j
0
c
i
j
0
0
-
c
i
j
0
c
i
j
)
(
x
i
y
i
x
j
y
j
)
=
(
b
x
i
b
y
i
b
x
j
b
y
j
)
,
where b
x
i and b
y
j are the components of the force exerted at node i in x and y directions, respectively, and x
i and y
i are the coordinates of node i in x and y directions, respectively. The spring systems work as a two-dimensional truss such that nodes and springs are points and elements of truss, respectively. We mention that the system of algebraic equations (22) can be assembled in its standard finite element concept to yield global force vector for whole system of springs as follows:
(23)
C
X
=
B
,
where C is the stiffness matrix of the system calculated by assembling the stiffness matrices of all the springs defined in the system and B represents the vector of nodal forces. In the equilibrium condition, the vector B of assembled spring forces should be equal to zero. This requirement leads to the following system of algebraic equation which should be solved for the unknown vector of nodal position X, that is, solve
(24)
C
X
=
0
.
It is obvious that the equation system defined in (24) is singular before any boundary conditions are considered. The boundary conditions used here for solving this system of equations are defined by the requirement that the boundary nodes should not be allowed to move perpendicularly to the boundaries. In other words, boundary nodes only can be displaced along the boundaries which they have been placed on. Mathematical representation of these boundary conditions can be defined as
(25)
Δ
X
i
T
n
i
=
(
Δ
x
i
Δ
y
i
)
(
n
x
i
n
y
i
)
=
(
(
x
i
-
x
-
ℓ
)
(
y
i
-
y
-
ℓ
)
)
(
n
x
i
n
y
i
)
=
0
,
where x
-
ℓ and y
-
ℓ are the initial coordinates of boundary node i; x
i and y
i represent the displaced final position of node i and n is outward unit vector normal to the boundary at node i. This condition guarantees that the nodes initially located at the intersection of two boundary lines must remain on its initial position. Boundary conditions for a simple net of springs are shown in Figure 3.Figure 3
Boundary conditions of spring system.Upon solving the system of equations with appropriate boundary conditions, the refined position of nodes is obtained leading to substantial reduction of the local and global error of the numerical solution in the subsequent analysis. The efficiency and effectiveness of the proposed adaptive refinement technique is verified in the next section by its application to benchmark test examples in plane elasticity.


#### 3.2.2. Node Enrichment

With the meshless methods, the enrichment strategy only requires that the locations of new nodes to be added are determined without requiring to define the connectivity of the resulting configuration. Different methods can be thought to be defining the location of the nodes to be added to the current nodal configuration. Here the new nodes are added in the neighborhood of existing nodal points defined by a Voronoi diagram.Once the Voronoi cells are defined, the vertices of Voronoi cells corresponding to the nodes with higher error than the average error over the domain are considered as the new nodes to be added to the current nodal configuration. The method is schematically illustrated in Figure4. The value of average error e
avg over the domain is obtained by
(26)
e
avg
=
∑
ℓ
=
1
M
e
ℓ
M
.
The above node enrichment adaptive refinement procedure based on error estimation is schematically shown in Figure 5.Figure 4
The node enrichment using Voronoi diagram (∙ initial nodes, ▲ new added nodes).Figure 5
The node enrichment procedure in MDLSM method.


## 3.1. Error Estimator

In numerical methods, a problem is solved by discretization of the problem domain into the subdomains, hence, the governing equations only apply into these subdomains, so numerical methods always come with discretization error. Discretization error is one of the most important challenges in the numerical methods. Discretization error is theoretically decreased by refining the discretization domain but perfunctory refinement imposes the heavy computational cost without supplying the expected accuracy. Adaptive refinement methods mean balances between refinement procedure and its computational cost. These methods only refine locally the regions of the domain which has higher error. Adaptive procedure has two main parts: error estimation and adaptive refinement. Any success adaptive refinement needs a reliable error estimation procedure. Real error distribution can not be practicably used because the exact solution is not available for any practical problems. Several methods are used for error estimation with different numerical methods and these methods are categorized into two classes, namely, the residual based method [34] and recovery based method [35]. In residual based method, the residuals of differential equation and its boundaries are used as a criterion of error. The gradient of the solutions is used in recovery based method as the error criterion.In this paper, the relative least square residuals functional for each node is defined as follows:(18)
e
=
I
(
U
t
T
U
t
)
,
where I is the least square residuals functional in (15) and U
t is unknowns obtained from the main solution. It is noticed that most of the computations of the least square residuals functional can be obtained from the main solution of the MDLSM method.


## 3.2. Adaptive Refinement Procedures

### 3.2.1. Node Moving

Mesh movement strategy can be easily and efficiently used with meshless methods since no element distortion is associated with the method. It should be noted that the mesh movement technique can be used in conjunction with the MDLSM method to adaptively adjust nodal points to improve the quality of the solution obtained with a prespecified number of nodal points. Here, a nodal refinement procedure is used that is called node moving adaptive refinement approach. When a node refinement is required, springs of prescribed stiffness are placed between each pair of nodes belonging to the same subdomain and the nodes are then moved until the spring system is in equilibrium.In the node moving procedure, first, all nodes are connected with springs in which the neighbor nodes are defined using Voronoi diagram [36] (see Figure 1). Voronoi diagram is defined as
(19)
T
i
=
{
x
∈
ℝ
:
d
(
x
,
x
i
)
<
d
(
x
,
x
j
)
,
j
≠
i
}
,
where d
(
x
,
x
i
) is Euclidean distance between x and x
i. The above equation means that neighbor nodes to node i are the closest nodes to the node i rather than other nodes.Figure 1
Voronoi diagram of neighboring nodes and spring connections.Spring forces are defined as(20)
b
i
j
=
c
i
j
(
x
i
-
x
j
)
,
where c
i
j is stiffness of spring between i, j and x
i, x
j are coordinates of i, j in equilibrium, respectively. The free body diagram is shown in Figure 2.Figure 2
Free body diagram between two nodesi and j.Spring stiffness is defined as a function of errors between two pointsi, j as follows:
(21)
c
i
j
=
(
e
i
+
e
j
)
d
i
j
,
where e
i and e
j are the values of the error estimators obtained from (18) at nodes i and j, respectively, and d
i
j is distance between these two connected nodes. In matrix form we have
(22)
(
c
i
j
0
-
c
i
j
0
0
c
i
j
0
-
c
i
j
-
c
i
j
0
c
i
j
0
0
-
c
i
j
0
c
i
j
)
(
x
i
y
i
x
j
y
j
)
=
(
b
x
i
b
y
i
b
x
j
b
y
j
)
,
where b
x
i and b
y
j are the components of the force exerted at node i in x and y directions, respectively, and x
i and y
i are the coordinates of node i in x and y directions, respectively. The spring systems work as a two-dimensional truss such that nodes and springs are points and elements of truss, respectively. We mention that the system of algebraic equations (22) can be assembled in its standard finite element concept to yield global force vector for whole system of springs as follows:
(23)
C
X
=
B
,
where C is the stiffness matrix of the system calculated by assembling the stiffness matrices of all the springs defined in the system and B represents the vector of nodal forces. In the equilibrium condition, the vector B of assembled spring forces should be equal to zero. This requirement leads to the following system of algebraic equation which should be solved for the unknown vector of nodal position X, that is, solve
(24)
C
X
=
0
.
It is obvious that the equation system defined in (24) is singular before any boundary conditions are considered. The boundary conditions used here for solving this system of equations are defined by the requirement that the boundary nodes should not be allowed to move perpendicularly to the boundaries. In other words, boundary nodes only can be displaced along the boundaries which they have been placed on. Mathematical representation of these boundary conditions can be defined as
(25)
Δ
X
i
T
n
i
=
(
Δ
x
i
Δ
y
i
)
(
n
x
i
n
y
i
)
=
(
(
x
i
-
x
-
ℓ
)
(
y
i
-
y
-
ℓ
)
)
(
n
x
i
n
y
i
)
=
0
,
where x
-
ℓ and y
-
ℓ are the initial coordinates of boundary node i; x
i and y
i represent the displaced final position of node i and n is outward unit vector normal to the boundary at node i. This condition guarantees that the nodes initially located at the intersection of two boundary lines must remain on its initial position. Boundary conditions for a simple net of springs are shown in Figure 3.Figure 3
Boundary conditions of spring system.Upon solving the system of equations with appropriate boundary conditions, the refined position of nodes is obtained leading to substantial reduction of the local and global error of the numerical solution in the subsequent analysis. The efficiency and effectiveness of the proposed adaptive refinement technique is verified in the next section by its application to benchmark test examples in plane elasticity.


### 3.2.2. Node Enrichment

With the meshless methods, the enrichment strategy only requires that the locations of new nodes to be added are determined without requiring to define the connectivity of the resulting configuration. Different methods can be thought to be defining the location of the nodes to be added to the current nodal configuration. Here the new nodes are added in the neighborhood of existing nodal points defined by a Voronoi diagram.Once the Voronoi cells are defined, the vertices of Voronoi cells corresponding to the nodes with higher error than the average error over the domain are considered as the new nodes to be added to the current nodal configuration. The method is schematically illustrated in Figure4. The value of average error e
avg over the domain is obtained by
(26)
e
avg
=
∑
ℓ
=
1
M
e
ℓ
M
.
The above node enrichment adaptive refinement procedure based on error estimation is schematically shown in Figure 5.Figure 4
The node enrichment using Voronoi diagram (∙ initial nodes, ▲ new added nodes).Figure 5
The node enrichment procedure in MDLSM method.


## 3.2.1. Node Moving

Mesh movement strategy can be easily and efficiently used with meshless methods since no element distortion is associated with the method. It should be noted that the mesh movement technique can be used in conjunction with the MDLSM method to adaptively adjust nodal points to improve the quality of the solution obtained with a prespecified number of nodal points. Here, a nodal refinement procedure is used that is called node moving adaptive refinement approach. When a node refinement is required, springs of prescribed stiffness are placed between each pair of nodes belonging to the same subdomain and the nodes are then moved until the spring system is in equilibrium.In the node moving procedure, first, all nodes are connected with springs in which the neighbor nodes are defined using Voronoi diagram [36] (see Figure 1). Voronoi diagram is defined as
(19)
T
i
=
{
x
∈
ℝ
:
d
(
x
,
x
i
)
<
d
(
x
,
x
j
)
,
j
≠
i
}
,
where d
(
x
,
x
i
) is Euclidean distance between x and x
i. The above equation means that neighbor nodes to node i are the closest nodes to the node i rather than other nodes.Figure 1
Voronoi diagram of neighboring nodes and spring connections.Spring forces are defined as(20)
b
i
j
=
c
i
j
(
x
i
-
x
j
)
,
where c
i
j is stiffness of spring between i, j and x
i, x
j are coordinates of i, j in equilibrium, respectively. The free body diagram is shown in Figure 2.Figure 2
Free body diagram between two nodesi and j.Spring stiffness is defined as a function of errors between two pointsi, j as follows:
(21)
c
i
j
=
(
e
i
+
e
j
)
d
i
j
,
where e
i and e
j are the values of the error estimators obtained from (18) at nodes i and j, respectively, and d
i
j is distance between these two connected nodes. In matrix form we have
(22)
(
c
i
j
0
-
c
i
j
0
0
c
i
j
0
-
c
i
j
-
c
i
j
0
c
i
j
0
0
-
c
i
j
0
c
i
j
)
(
x
i
y
i
x
j
y
j
)
=
(
b
x
i
b
y
i
b
x
j
b
y
j
)
,
where b
x
i and b
y
j are the components of the force exerted at node i in x and y directions, respectively, and x
i and y
i are the coordinates of node i in x and y directions, respectively. The spring systems work as a two-dimensional truss such that nodes and springs are points and elements of truss, respectively. We mention that the system of algebraic equations (22) can be assembled in its standard finite element concept to yield global force vector for whole system of springs as follows:
(23)
C
X
=
B
,
where C is the stiffness matrix of the system calculated by assembling the stiffness matrices of all the springs defined in the system and B represents the vector of nodal forces. In the equilibrium condition, the vector B of assembled spring forces should be equal to zero. This requirement leads to the following system of algebraic equation which should be solved for the unknown vector of nodal position X, that is, solve
(24)
C
X
=
0
.
It is obvious that the equation system defined in (24) is singular before any boundary conditions are considered. The boundary conditions used here for solving this system of equations are defined by the requirement that the boundary nodes should not be allowed to move perpendicularly to the boundaries. In other words, boundary nodes only can be displaced along the boundaries which they have been placed on. Mathematical representation of these boundary conditions can be defined as
(25)
Δ
X
i
T
n
i
=
(
Δ
x
i
Δ
y
i
)
(
n
x
i
n
y
i
)
=
(
(
x
i
-
x
-
ℓ
)
(
y
i
-
y
-
ℓ
)
)
(
n
x
i
n
y
i
)
=
0
,
where x
-
ℓ and y
-
ℓ are the initial coordinates of boundary node i; x
i and y
i represent the displaced final position of node i and n is outward unit vector normal to the boundary at node i. This condition guarantees that the nodes initially located at the intersection of two boundary lines must remain on its initial position. Boundary conditions for a simple net of springs are shown in Figure 3.Figure 3
Boundary conditions of spring system.Upon solving the system of equations with appropriate boundary conditions, the refined position of nodes is obtained leading to substantial reduction of the local and global error of the numerical solution in the subsequent analysis. The efficiency and effectiveness of the proposed adaptive refinement technique is verified in the next section by its application to benchmark test examples in plane elasticity.


## 3.2.2. Node Enrichment

With the meshless methods, the enrichment strategy only requires that the locations of new nodes to be added are determined without requiring to define the connectivity of the resulting configuration. Different methods can be thought to be defining the location of the nodes to be added to the current nodal configuration. Here the new nodes are added in the neighborhood of existing nodal points defined by a Voronoi diagram.Once the Voronoi cells are defined, the vertices of Voronoi cells corresponding to the nodes with higher error than the average error over the domain are considered as the new nodes to be added to the current nodal configuration. The method is schematically illustrated in Figure4. The value of average error e
avg over the domain is obtained by
(26)
e
avg
=
∑
ℓ
=
1
M
e
ℓ
M
.
The above node enrichment adaptive refinement procedure based on error estimation is schematically shown in Figure 5.Figure 4
The node enrichment using Voronoi diagram (∙ initial nodes, ▲ new added nodes).Figure 5
The node enrichment procedure in MDLSM method.


## 4. Numerical Experiments

In this section, we are solving the benchmark examples by using the proposed node moving and node enrichment refinement procedures for the MDLSM method and comparing the results with the refinement procedure in the irreducible DLSM method and the exact analytical solutions or finite element results with very fine mesh.Example 1 (an infinite plate with a circular hole).
In the first example, consider the case of an infinite plate with a circular hole subjected to a uniaxial tractionP at infinity, as shown in Figure 6. Due to symmetry, only the upper right square quadrant of the plate is modeled. The edge length of the square is 5
a, where a is the radius of the circular hole. This example is chosen because the exact analytical solution is available from Timoshenko and Goodier [37]. The solutions for the displacements and the stresses under a unit uniaxial stress along the x axis are given as follows:Figure 6
An infinite plate with a circular hole under a uniaxial loadP (Example 1).(27)
σ
x
=
t
(
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
=
-
t
(
a
2
r
2
(
1
2
sin
(
2
θ
)
+
sin
(
4
θ
)
)
-
3
a
4
2
r
4
sin
(
4
θ
)
)
,
u
r
=
t
4
G
(
r
(
κ
-
1
2
+
cos
⁡
(
2
θ
)
)
+
a
2
r
(
1
+
(
1
+
κ
)
cos
⁡
(
2
θ
)
)
-
a
4
r
3
cos
⁡
(
2
θ
)
)
,
u
θ
=
t
4
G
(
(
1
-
κ
)
a
2
r
-
r
-
a
4
r
3
)
sin
(
2
θ
)
,whichG is the shear modulus and κ
=
(
3
-
ν
)
/
(
1
+
ν
) where ν is the Poissons ratio. In this example, the constant values are a
=
1, t
=
1, and E
=
1000 and ν
=
0.3.For the node moving procedure, initial and refined nodal configurations with 122, 229, and 305 nodes are shown in Figure7. The nodal points are refined with respect to the proposed error estimator based on the least square residuals functional. In Figure 8(a) numerical results of normal stress σ
x on the left edge is compared with the exact analytical solution. Figures 8(b) and 8(c) compare the convergence curve and the convergence rate of the MDLSM method for the node moving adaptive refinement strategy and plotted using the error norms in Table 1.Table 1
Comparison of the error norms for initial and refined (using node moving) nodal configurations (Example1).
Number of nodes
Norm of error for initial nodes
Norm of error for refined nodes
122
0.6684
0.1528
229
0.3550
0.0331
305
0.2123
0.0206Initial and refined (node moving) nodal configurations (Example1).
(a)
122 initial nodes
(b)
122 refined nodes
(c)
229 initial nodes
(d)
229 refined nodes
(e)
305 initial nodes
(f)
305 refined nodes(a) Normal stressσ
x at x
=
0. (b) Convergence curve of node moving procedure. (c) Convergence rate of node moving procedure.
(a)
(b)
(c)Also, the node enrichment adaptive refinement strategy is applied to the MDLSM method for obtaining accurate results. As shown in Figure9, first the problem domain is discretized by using 122 initial nodal distributions, and then the nodal points are enriched in three steps (182, 288, and 521 nodes) on the region determined by the error indicator. Figure 10(a) compares the numerical results of normal stress σ
x on the left edge with exact analytical solution. The error norms in Table 2 are used for plotting and comparison of the convergence curve and the convergence rate of the node enrichment procedure for the MDLSM method (see Figures 12(b) and 12(c)). The results clearly show that the node moving and multistage node enrichment adaptive refinement strategies in the MDLSM method are more efficient compared to the refinement procedures in the DLSM method. Figures 8(b) and 8(c), 10(b), and 10(c) indicate that by using the node moving and the node enrichment adaptive refinement procedures, the convergence rate of the MDLSM method is increased.Table 2
Error norms in initial nodal configuration and in different steps of multistage node enrichment procedure (Example1).
Number of nodes
Norm of the error
122
0.6684
182
0.1457
288
0.0571
521
0.0190Initial and refined (node enrichment) nodal configurations (Example1).
(a)
Initial nodes = 122
(b)
First step enrichment nodes = 182
(c)
Second step enrichment nodes = 288
(d)
Third step enrichment nodes = 521(a) Normal stressσ
x
x at x
=
0. (b) Convergence curve for node enrichment procedure. (c) Convergence rate for node enrichment procedure.
(a)
(b)
(c)Example 2 (a cylinder subjected to an internal pressure).
As a second elastostatic benchmark example a cylinder subjected to an internal pressure is considered. Due to the symmetry, only a quarter of the cylinder is modeled; see Figure11. The boundary conditions are illustrated in Figure 11. The exact analytical solution of this problem is
(28)
σ
r
=
a
2
P
b
2
-
a
2
(
1
-
b
2
r
2
)
,
σ
θ
=
a
2
P
b
2
-
a
2
(
1
+
b
2
r
2
)
,
where the constant values are a
=
1, b
=
5, P
=
1, ν
=
0.3, and E
=
1
0
7. As shown in Figure 12, four types of nodal distributions with 106, 214, 359, and 543 nodes are distributed to solve and refine the nodes in the problem domain.Figure 11
A cylinder subjected to an internal pressure and its boundary conditions.Initial and refined (node moving) nodal configurations (Example2).
(a)
106 initial nodes
(b)
106 refined nodes
(c)
214 initial nodes
(d)
214 refined nodes
(e)
359 initial nodes
(f)
359 refined nodes
(g)
543 initial nodes
(h)
543 refined nodesFigure13(a) compares the normal stress σ
x at x
=
0 for initial and refined nodal configurations with 543 nodes. It is clear that the result of refined nodal configuration is more similar to the exact analytical solution than initial nodal configuration. Table 3 compares the error norms of the node moving procedure for the initial and refined nodal configurations and Figures 13(b) and 13(c) compare the convergence curve and the convergence rate, respectively.Table 3
Comparison of the error norms for initial and refined (using node moving) nodal configurations (Example2).
Number of nodes
Norm of error for initial nodes
Norm of error for refined nodes
106
0.3843
0.1567
214
0.1247
0.0388
359
0.0691
0.0153
543
0.0354
0.0098(a) Normal stressσ
x at x
=
0. (b) Convergence curve of node moving procedure. (c) Convergence rate of node moving procedure.
(a)
(b)
(c)In Figure14, the initial and refined nodal distributions for the node enrichment strategy are shown. Figure 15(a) compares the normal stress σ
x at x
=
0 for initial and last steps refined with 543 nodal distributions. In Table 4 the error norm of the node enrichment refinement procedure is shown. Figures 15(b) and 15(c) compare the convergence curve and the convergence rate of the node enrichment adaptive refinement strategy.Table 4
Error norms in initial nodal configuration and in different steps of multistage node enrichment procedure (Example2).
Number of nodes
Norm of the error
106
0.3843
164
0.0808
276
0.0501
491
0.0191Initial and refined (node enrichment) nodal configurations (Example2).
(a)
Initial nodes = 106
(b)
First step enrichment nodes = 164
(c)
Second step enrichment nodes = 276
(d)
Third step enrichment nodes = 491(a) Normal stressσ
x at x
=
0. (b) Convergence curve for node enrichment procedure. (c) Convergence rate for node enrichment procedure.
(a)
(b)
(c)Example 3 (a reservoir fully filled with water).
In this example, consider that the wall of a reservoir fully filled with water is investigated. The geometry of the wall is irregular as given in Figure16. The material properties of the wall are given as Young’s modulus E
=
1
0
7 and Poissons ratio ν
=
0.3. The bottom of the wall is fixed and the curved edge of the wall is subjected to a hydrostatic pressure P
=
-
9800
(
H
-
y
)
MPa. Since the analytical solution of this problem is not available, a very fine mesh (with 59,400 linear triangular elements) FEM solution will be considered as our reference solution.Figure 16
A reservoir fully filled with water.As shown in Figure17, four types of nodal distribution with 84, 138, 218, and 299 nodes are distributed to solve and refine the points in the problem domain. Figure 18(a) compares the displacement in y-direction along the vertical edge for initial and refined configurations. Table 5 compares the error norms of the node moving procedures based on MDLSM method and Figures 18(a) and 18(b) compare the convergence curve and convergence rate, respectively. For node enrichment adaptive refinement procedure, the problem domain is discretized with initial 84 points and is refined in three steps near high gradient error norm and is solved with 120, 181, and 284 nodes, respectively (see Figure 19). Figure 20(a) compares the displacement in y-direction along the vertical edge for initial and refined configurations. In Table 6 the error norms of the node enrichment refinement procedure is used to plot the convergence curve and convergence rate in Figures 20(b) and 20(c).Table 5
Comparison of the error norms for initial and refined (using node moving) nodal configurations (Example3).
Number of nodes
Norm of error for initial nodes
Norm of error for refined nodes
84
0.0195
0.0190
138
0.0135
0.0116
218
0.0091
0.0059
299
0.0071
0.0047Table 6
Error norms in initial nodal configuration and in different steps of multistage node enrichment procedure (Example3).
Number of nodes
Norm of the error
84
0.0195
120
0.0147
181
0.0087
284
0.0060Initial and refined (node moving) nodal configurations (Example3).
(a)
84 initail nodes
(b)
84 refined nodes
(c)
138 initial nodes
(d)
138 refined nodes
(e)
218 initial nodes
(f)
218 refined nodes
(g)
299 initial nodes
(h)
299 refined nodes(a) Displacement iny-direction along the vertical edge. (b) Convergence curve of node moving procedure. (c) Convergence rate of node moving procedure.
(a)
(b)
(c)Initial and refined (node enrichment) nodal configurations (Example3).
(a)
Initial nodes = 84
(b)
First step enrichment nodes = 120
(c)
Second step enrichment nodes = 181
(d)
Third step enrichment nodes = 284(a) Displacement iny-direction along the vertical edge. (b) Convergence curve for node enrichment procedure. (c) Convergence rate for node enrichment procedure.
(a)
(b)
(c)


## 5. Conclusion

A mixed discrete least squares meshless method was extended for node moving and node enrichment adaptive refinements for efficient analysis of the planar elasticity problem. For the refinement procedures an error estimator based on least square residuals functional was formulated and used. Voronoi diagram was extended in the refinement procedures to find the neighbor nodes (node moving) and the position of the new nodes (node enrichment). For the moving node procedure, spring analogy was used to construct a system for computing the new place of each node after the refinement procedure. The efficiency and effectiveness of the proposed node moving and node enrichment adaptive refinement techniques in the MDLSM method by their application to the benchmark examples in the elasticity problems were verified. Results show that the proposed refinement methods are accurate and straightforward.



---

*Source: 721240-2014-02-16.xml*