An interpolating boundary element-free method for three-dimensional potential problems

Applied Mathematical Modelling 39 (2015) 3116–3134

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

An interpolating boundary element-free method for three-dimensional potential problems Xiaolin Li College of Mathematics Science, Chongqing Normal University, Chongqing 400047, PR China

a r t i c l e

i n f o

Article history: Received 12 October 2012 Received in revised form 21 August 2013 Accepted 23 October 2014 Available online 27 November 2014 Keywords: Meshless method Improved interpolating moving least-square method Boundary integral equations Interpolating boundary element-free method

a b s t r a c t This paper begins by discussing an improved interpolating moving least-square (IIMLS) method and the properties of its shape function. In the IIMLS method, the shape function is of delta function property and the weight function is nonsingular, so it overcomes the drawbacks in both the moving least-square approximation and the interpolating moving least-square method. Then combining boundary integral equations with the IIMLS method, an interpolating boundary element-free method (IBEFM) is developed for three-dimensional potential problems. In the IBEFM, only a nodal data structure on the boundary face of a domain is required. Unlike the boundary node method (BNM), the IBEFM is a direct meshless method in which the primary unknown quantities are real solutions of nodal variables, and boundary conditions can be imposed directly and easily, which leads to a greater computational precision. Besides, the number of both unknowns and system equations in the IBEFM is only half of that in the BNM, and thus the computing speed and efficiency are increased. Numerical examples on curve/surface fittings and potential problems indicate that the efficiency and convergence rate of the present methods is high. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction The finite element method (FEM) and the boundary element method (BEM) have been the dominant numerical computational methods in the field of computational science and engineering for several decades. Both methods depend on the generation of meshes, adapted or not. Mesh generation in some situations is still arduous, time consuming and fraught with pitfalls. To circumvent the problems associated with meshing, meshless (or meshfree) methods to obtain numerical solutions of boundary value problems without resorting to an element frame have been developed in the past two decades. The moving least-square (MLS) [1] is an approximation method, developed by Lancaster and Salkauskas, which generates continuous functions from a cluster of unorganized sampled point values based on the computation of a weighted least squares approximation. Li and Zhu [2,3] obtained error estimates of the MLS approximation in Sobolev spaces. Since the numerical approximations of the MLS start from scattered nodes instead of meshes, there have many meshless methods based on the MLS scheme. Among them are the diffuse element method, the element-free Galerkin (EFG) method [4], the h-p meshless method [5], the meshless local Petrov–Galerkin (MLPG) method [6], the moving least square reproducing kernel method [7,8], the reproducing kernel hierarchical partition of unity [9,10], and so on. These meshless methods followed the idea as the FEM, in which the problem domain is discretized.

E-mail address: [email protected] http://dx.doi.org/10.1016/j.apm.2014.10.071 0307-904X/Ó 2014 Elsevier Inc. All rights reserved.

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

3117

The MLS approximation has also been used in boundary integral equations (BIEs). Typical of them are the meshless local boundary integral equation (LBIE) method [6] and the boundary node method (BNM) [11–13]. As the BEM, these BIEs-based meshless methods have emerged as promising numerical techniques in scientific computing. Nonetheless, since the MLS approximations lack the delta function property of the usual FEM and BEM shape functions, boundary conditions in these meshless techniques are difficult to be implemented. The technique used in the BNM involves a new definition of the discrete norm used for the construction of the MLS approximations, which doubles the number of unknowns and system equations. Accordingly, the advantage of BIEs-based meshless methods is therefore eroded and discounted to a certain extent. via combining a variational form of BIEs and the MLS approximations, another technique is developed in the symmetric Galerkin BNM to impose boundary conditions [14–17]. However, in all these BIEs-based meshless methods, the basic unknown quantities are approximations of nodal variables, and therefore they are indirect meshless methods. Recently, Liew et al. [18] developed an improved MLS scheme that uses weighted orthogonal polynomials as basis functions. The improved MLS scheme has been introduced into BIEs to develop a direct meshless method, the boundary element-free method (BEFM) [18,19]. Because the improved MLS scheme still lacks the delta function property, boundary conditions in the BEFM are implemented with constraints [20,21]. To restore the delta function property of the MLS, Lancaster and Salkauskas further developed an interpolating moving least-square (IMLS) method that uses specific singular functions as weight functions [1]. Based on the IMLS scheme, an improved EFG has been proposed by Kaljevic and Saigal [22]. Besides, by revising the formulae of the IMLS method and combining it with BIEs, Ren et al. proposed an interpolating boundary element-free method (IBEFM) [20,21] and an interpolating element-free Galerkin (IEFG) method [23,24] for two-dimensional potential and elasticity problems. In these improved schemes, boundary conditions can be imposed as easily as in the FEM and the BEM. A drawback of the IMLS method is that the involved weight function is singular at nodes. The singularity complicates the computation of the inverse of the singular matrix and thus reduces the computing speed and efficiency. To avoid this shortcoming, Netuzhylov [25] introduced a perturbation technique into the IMLS method. This method requires a new parameter e, which is used for regularization of the weight function matrix. Although the parameter will affect the performance of this method, there is still no theoretical strategy to get appropriate values for this parameter. In order to overcome problems in both the MLS approximation and the IMLS method, Wang et al. [26,27] proposed an improved interpolating moving least-square (IIMLS) method, and an improved interpolating EFG method and an IBEFM for two-dimensional potential problems. The IIMLS method has the following features: (1) the shape function possesses the delta function property, so implementing boundary conditions in the IIMLS-based meshless method is much easier than that in the MLS-based meshless method; (2) the weight function used is nonsingular at any point, hence any weight function used in the MLS approximation can also be used in the IIMLS method; and (3) the number of unknown coefficients in the trial function of the IIMLS method is less than that of the MLS approximation, thus fewer nodes are required in the influence domain and a higher computational accuracy can be achieved in the IIMLS method. In this paper, we first discuss the IIMLS method and the properties of its shape function. Specially, the reproducing properties of the IIMLS shape function is proved theoretically and verified numerically in detail. The capability of data fittings of the IIMLS method for one- and two-dimensional functions is also studied. Then, a new implementation of the interpolating boundary element-free method (IBEFM) is developed for boundary-only analysis of three-dimensional potential problems. In our implementation, the IBEFM combines the IIMLS for construction of interpolation functions with BIEs for the governing partial differential equations. Besides, the IIMLS interpolants have been suitably truncated at corners and edges in order to avoid the discontinuity of the IIMLS shape function. Since potential problem is one of the most important problems in mathematics which has wide applications to a number of topics relevant to mathematical physics and engineering, the current IBEFM is formulated for three-dimensional potential problems [28]. In the IBEFM, the boundary of the problem domain is discretized via properly scattered nodes. Since the IIMLS shape functions satisfy delta function properties, the IBEFM overcomes the shortcomings of the BNM. The enforcement of boundary conditions in the IBEFM is as easy as in the conventional BEM. Therefore, the number of both system equations and unknowns in the IBEFM is only half of that in the BNM. The rest of this paper is outlined as follows. Section 2 presents the IIMLS method and the properties of its shape function. Curve and surface fittings as examples are also studied in this section to show the capacity and accuracy of the IIMLS method. In Section 3, a detailed numerical implementation of the IBEFM is described. A comparison between the IBEFM, the BNM and the BEM is also included in this section. Then, numerical results are given in Section 4. Section 5 contains conclusions. 2. The improved interpolating moving least-squares (IIMLS) method 2.1. Notations Let D be a bounded domain in Rd (d ¼ 2; 3Þ. A general point of Rd is denoted as x ¼ ðx1 ; x2 ; . . . ; xd Þ or y ¼ ðy1 ; y2 ; . . . ; yd Þ. For any x 2 D, assume that the influence domain of x is a d-dimensional ball RðxÞ with radius rðxÞ, i.e.,

RðxÞ ¼ fy : kx  yk 6 r ðxÞg; where kk may denote any norm in Rd . For simplicity, kk is the Euclidean norm. Besides, given x 2 D, assume that there have n nodes xi such that x 2 Ri ; i ¼ 1; 2; . . . ; n, where Ri ¼ Rðxi Þ is the influence domain of xi .

3118

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

Moreover, let wi ðxÞ :¼ wðx  xi Þ be a weight function that satisfies:

x 2 Ri () wi ðxÞ > 0 and wi ðxÞ – 0;

ð1Þ

for all x 2 D and i ¼ 1; 2; . . . ; n. 2.2. The IIMLS procedure Let fpi ðxÞgm i¼0 denote a given set of basis functions with p0 ðxÞ  1. The commonly used are polynomial basis functions. In this study, we take quadratic polynomial basis. Before generating shape functions with delta function property, we first construct the following basis functions,

Þ ¼ pi ðx Þ  g i ðx; x

n X     q x; xj pi xj ;

i ¼ 0; 1; . . . ; m;

ð2Þ

j¼1

   2 RðxÞ can either be the evaluation point x or a nodal point xj . The function q x; xj is chosen such that where x

  q xi ; xj ¼ dij ;

ð3Þ

and n X   q x; xj ¼ 1;

8x 2 D:

ð4Þ

j¼1

Similar to the function used in Ref. [26], this function is taken as





q x; xj ¼

8 < Pnnðx;xj Þ :

k¼1

nðx;xk Þ

0;

; x 2 Rj ;

1 6 j 6 n;

ð5Þ

x R Rj ;

with

Qn s   k¼1;k–j kx  xk k   : n x; xj ¼ Qn xj  xk s k¼1;k–j Here, s > 0 and kk may denote any norm in Rd , but for reasons of differentiability it is natural to employ Euclidean norm. Þ  0 for any x 2 D. In Eq. (2) letting i ¼ 0, then noting that p0 ðxÞ  1 and using Eq. (4) we obtain g 0 ðx; x Þ as ~ ðx; x To obtain the IIMLS scheme for a given function v ðxÞ, we first define a new function v

v~ ðx; xÞ ¼ v ðfxÞ 

n X     q x; xj v xj :

s ¼ 2 and the

ð6Þ

j¼1

 Þ of the function v Þ at a point x is defined in RðxÞ by ~ h ðx; x ~ ðx; x Then, the local approximation v

v~ h ðx; xÞ ¼

m X Þaj ðxÞ: g j ðx; x

ð7Þ

j¼1

where aj ðxÞ are unknown coefficients and are obtained at any point x by minimizing the following weighted discrete L2 norm,

( " #)2 n n m n X X X X 2 J¼ wi ðxÞ½v~ h ðx; xi Þ  v~ ðx; xi Þ ¼ wi ðxÞ g j ðx; xi Þaj ðxÞ  v ðxi Þ  qðx; xj Þv ðxj Þ i¼1

i¼1

j¼1

j¼1

( )2 n m n X X X   ¼ wi ðxÞ g j ðx; xi Þaj ðxÞ  dij  qðx; xj Þ v ðxj Þ : i¼1

j¼1

j¼1

To find the coefficients aj ðxÞ, we obtain the extremum of J by

( ) n m n X X X      @J ¼ 2 wi ðxÞ g j ðx; xi Þaj ðxÞ  dij  q x; xj v xj g k ðx; xi Þ ¼ 0; @ak i¼1 j¼1 j¼1 Hence

k ¼ 1; 2; . . . ; m:

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

" # m n n n X X X X   wi ðxÞg j ðx; xi Þg k ðx; xi Þ aj ðxÞ ¼ wi ðxÞg k ðx; xi Þ dij  qðx; xj Þ v ðxj Þ; j¼1

i¼1

i¼1

3119

k ¼ 1; 2; . . . ; m:

j¼1

It is expressed in matrix form as follows:

AðxÞaðxÞ ¼ BðxÞv ;

ð8Þ

in which

AðxÞ ¼

n X wi ðxÞgðx; xi ÞgT ðx; xi Þ;

ð9Þ

i¼1

aðxÞ ¼ ½a1 ðxÞ; a2 ðxÞ; . . . ; am ðxÞT ; BðxÞ ¼

n X wi ðxÞgðx; xi Þ½ei  qðxÞ;

ð10Þ

i¼1

v ¼ ½v ðx1 Þ; v ðx2 Þ; . . . ; v ðxn ÞT : Here, ei is a n unit row vector with ith component being 1, and T gðx; xi Þ ¼ ½g 1 ðx; xi Þ; g 2 ðx; xi Þ; . . . ; g m ðx; xi Þ ;

ð11Þ

qðxÞ ¼ ½qðx; x1 Þ; qðx; x2 Þ; . . . ; qðx; xn Þ: Unique solution of Eq. (8) is obtained if the inverse of matrix AðxÞ exists,

aðxÞ ¼ A1 ðxÞBðxÞv ;

ð12Þ

Substituting Eq. (12) into Eq. (7) yields the local approximation

v~ h ðx; xÞ ¼ gT ðx; xÞA1 ðxÞBðxÞv :

ð13Þ

 ¼ x, then from Eqs. (6) and (13), the global approximation of Setting x n n X X v h ðxÞ ¼ wi ðxÞv ðxi Þ þ qðx; xi Þv ðxi Þ: i¼1

v ðxÞ can be represented as ð14Þ

i¼1

where

wi ðxÞ ¼

m X

h i g j ðxÞ A1 ðxÞBðxÞ ;

ð15Þ

ji

j¼1

Here, g j ðxÞ is defined by Eq. (2) as

g j ðxÞ ¼ pj ðxÞ 

n X qðx; xk Þpj ðxk Þ;

j ¼ 1; 2; . . . ; m:

ð16Þ

k¼1

We define the IIMLS shape functions as

Ui ðxÞ ¼ wi ðxÞ þ qðx; xi Þ;

i ¼ 1; 2; . . . ; n;

then the global approximation of

v h ðxÞ ¼

n X

Ui ðxÞv ðxi Þ:

ð17Þ

v ðxÞ is ð18Þ

i¼1

Remark 2.1. For the construction of the IIMLS shape functions, the m  m moment matrix AðxÞ needs to be invertible. In the MLS approximations [2,3,5], a necessary condition for the invertibility of the moment matrix is that there are at least m þ 1 nodes covered in RðxÞ, where m þ 1 is the number of basis functions employed in the approximation. Similarly, a necessary condition in the IIMLS approximations is that there are at least m nodes covered in RðxÞ. Besides, Eq. (7) indicates that the number of the unknown coefficients aj ðxÞ in the IIMLS method is fewer than that in the MLS approximation. As a result, fewer nodes can be chosen in the meshless method formed with the IIMLS method compared to that formed with the MLS approximation.

Remark 2.2. In the IMLS method developed by Lancaster and Salkauskas [1], only singular weight function can be used. In contrast, any weight function used in the MLS approximation can also be used in the current IIMLS method. A variety of

3120

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

weight functions have been extensively investigated in the past for meshless methods. It has been shown in Refs. [12,13] that the Gaussian weight function yields excellent results. Hence, the numerical results in this work have been obtained using only the following Gaussian weight function:

( wi ðsÞ ¼

2

e½di =ð^cqi Þ ; 0 6 di 6 qi ; di P qi ;

0;

where di is the distance between an evaluation point x and the node xi ; ci is a constant controlling the shape of the weight function, qi is the size of the compact support of the weight function. Remark 2.3. The point interpolation method (PIM) [29] is another meshless interpolation technique. In the PIM, the polynomial basis and the radial basis are used for generating meshless shape functions with delta function properties. Some BIEsbased meshless methods using the PIM shape functions, such as the boundary PIM [30,31] and the hybrid radial BNM [32,33] have been developed. Boundary conditions in these meshless methods can also be imposed directly and easily. 2.3. Properties of IIMLS shape functions Proposition 2.1. The IIMLS shape function Ui ðxÞ has compact support Ri ; i ¼ 1; 2; . . . ; n. Proof. Eqs. (1), (9) and (10) imply that wi ðxÞ given in Eq. (15) has compact support Ri , while Eq. (5) implies that qðx; xi Þ also has compact support Ri . Thus, the proof is followed using Eq. (17). h Proposition 2.2. The IIMLS shape function Ui ðxÞ is of delta function properties, namely,



 

Ui xj ¼ dij ¼

1; i ¼ j; 0; i – j;

i; j ¼ 1; 2; . . . ; n:

ð19Þ

Proof. According to Eq. (3), it follows from Eq. (16) that n     X       g k xj ¼ pk xj  q xj ; xl pk ðxl Þ ¼ pk xj  pk xj ¼ 0; l¼1

for all k ¼ 2; . . . ; m and j ¼ 1; 2; . . . ; n. Then from Eq. (15) we have m   X  h    i wi xj ¼ g k xj A1 xj B xj ¼ 0: ki

k¼1

Finally, the proof is completed by using Eq. (17) and recalling again Eq. (3). h Fig. 1 shows a typical shape function for the MLS method and the IIMLS method in one-dimensional domain, while Fig. 2 shows a shape function in two-dimensional domain. Both figures indicate that the MLS shape functions do not satisfy delta function properties, while the IIMLS shape functions possess delta function properties. Proposition 2.3. The IIMLS shape function Ui ðxÞ is of reproducing properties as n X

Ui ðxÞpl ðxi Þ ¼ pl ðxÞ;

l ¼ 0; 1; . . . ; m:

ð20Þ

i¼1

In particular, Ui ðxÞ is of unity partition as n X

Ui ðxÞ ¼ 1:

ð21Þ

i¼1

Proof. From Eq. (17), n X

n n X X wi ðxÞpl ðxi Þ þ qðx; xi Þpl ðxi Þ:

i¼1

i¼1

Ui ðxÞpl ðxi Þ ¼

Using Eq. (15) yields

i¼1

ð22Þ

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

Fig. 1. Shape function of (a) the MLS method and (b) the IIMLS method for one-dimensional domain.

Fig. 2. Shape function of (a) the MLS method and (b) the IIMLS method for two-dimensional domain.

! n n X m n X m m h h i i X X X X wi ðxÞpl ðxi Þ ¼ g j ðxÞ A1 ðxÞBðxÞ pl ðxi Þ ¼ g j ðxÞ A1 ðxÞ ½BðxÞki pl ðxi Þ i¼1

ji

i¼1 j¼1

¼

m X m X n h X

A1 ðxÞ

j¼1 k¼1 i¼1

i jk

i¼1 j¼1

k¼1

½BðxÞki pl ðxi Þg j ðxÞ:

Eq. (10) shows

½BðxÞki ¼ wi ðxÞg k ðx; xi Þ  qðx; xi Þ

n X wa ðxÞg k ðx; xa Þ;

a¼1

thus n X wi ðxÞpl ðxi Þ ¼ M 1  M2 i¼1

where

M1 ¼

m X m X n h X

A1 ðxÞ

i jk

j¼1 k¼1 i¼1

M2 ¼

m X m X n h X

A1 ðxÞ

j¼1 k¼1 i¼1

¼

m X m h X j¼1 k¼1

Eq. (2) implies

i jk

wi ðxÞg k ðx; xi Þpl ðxi Þg j ðxÞ; n X qðx; xi Þ wa ðxÞg k ðx; xa Þpl ðxi Þg j ðxÞ

a¼1

( ) n n i X X 1 A ðxÞ wa ðxÞg k ðx; xa Þ qðx; xi Þpl ðxi Þ g j ðxÞ: jk

a¼1

i¼1

jk

3121

3122

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134 n X qðx; xi Þpl ðxi Þ ¼ pl ðxa Þ  g l ðx; xa Þ; i¼1

hence according to Eq. (9),

M2 ¼

m X m h X

n m X m h i X i X A ðxÞ wa ðxÞg k ðx; xa Þpl ðxa Þg j ðxÞ  A1 ðxÞ 1

jk

j¼1 k¼1

¼ M1 

m X m h X

a¼1

jk

a¼1

i A1 ðxÞ ½AðxÞkl g j ðxÞ: jk

j¼1 k¼1

Pm h

Then with the relation

j¼1 k¼1

( ) n X wa ðxÞg k ðx; xa Þg l ðx; xa Þ g j ðxÞ

k¼1

A1 ðxÞ

i jk

½AðxÞkl ¼ djl , we have

M2 ¼ M1  g l ðxÞ; so n X wi ðxÞpl ðxi Þ ¼ g l ðxÞ: i¼1

Finally, substituting the above equation into Eq. (22) and using Eq. (16) we can prove Eq. (20) as n X

n X qðx; xi Þpl ðxi Þ ¼ pl ðxÞ:

i¼1

i¼1

Ui ðxÞpl ðxi Þ ¼ g l ðxÞ þ

The proof of Eq. (21) follows immediately from the fact p0 ðxÞ  1.

h

2.4. Curve and surface fitting In this subsection, curve and surface fitting capability of the present IIMLS method is computed through one- and twodimensional functions. Let N denote the total number of nodes in the whole computational domain D. The analytical function v ðxÞ sampled at these nodes xi to determine the data v ðxi Þ; i ¼ 1; 2; . . . ; N. For the purpose of error estimation and convergence studies, a norm error is defined as [12,13]:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u Ns

2 X 100 u t1 eðv Þ ¼ v ðiÞ  v ðhiÞ percent; jv jmax Ns i¼1 ðiÞ

where N s is the number of sample points; v ðiÞ and v h refer to the values at the ith sample point obtained by the analytical solution and numerical method, respectively; jv jmax is the maximum value of v over N s points. 2.4.1. Patch test: reproducibility To confirm the reproducing properties of the IIMLS method given in Proposition 2.3, the standard patch test is conducted. Following two one-dimensional functions are test in the domain ½1; 1:

v 1 ðx1 Þ ¼ x1 ;

v 2 ðx1 Þ ¼ x21 :

Following two two-dimensional functions are test in the domain ½1; 1  ½1; 1:

v 3 ðx1 ; x2 Þ ¼ x1 þ x2 ;

v 4 ðx1 ; x2 Þ ¼ x21 þ x1 x2 þ x22 :

In this analysis, the quadratic polynomial basis is used. Besides, the one-dimensional functions are measured at 8 equidistant points, while the two-dimensional functions are measured at 16 equidistant points. The norm error for v 1 ðx1 Þ; v 2 ðx1 Þ; v 3 ðx1 ; x2 Þ and v 4 ðx1 ; x2 Þ is 1:208  1016 ; 1:488  1016 ; 1:707  1016 and 9:590  1017 , respectively. These numerical results imply that the analytical values are reproduced to machine precision. Therefore, the present IIMLS method passes the standard patch test exactly, which confirms Proposition 2.3. 2.4.2. Curve fitting Consider the function [1]

v ðx1 Þ ¼ ex

1

;

x1 2 ½5; 5:

Fig. 3 gives the log–log plot of errors with respect to the number of nodes. It can be seen that, as the number of nodes increases, the error becomes smaller. The experimental convergence rate is very high and approximate to 3. In Fig. 3, the results of the MLS method are also presented for comparison. It can be found that the errors from the IIMLS method are less than that from the MLS method. Besides, when a large number of nodes are used, the error of the MLS method does not decrease along with the increase of the nodes.

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

3123

Fig. 3. Convergence for curve fitting.

Fig. 4. Convergence for surface fitting.

2.4.3. Surface fitting Consider the function [34]

v ð x1 ; x2 Þ ¼



2 x21  x22 ;

x1 ; x2 2 ½1; 1:

The errors of the MLS method and the IIMLS method are shown in Fig. 4. The errors of both methods decrease along with the increase of the nodes. Again, high rates of convergence are achievable, and the present IIMLS method is more accurate. 2.4.4. Non-smooth function fitting In order to verify the power of the interpolation method and the significance of the strong locality, we choose the following specific function [34] on a rectangular domain D ¼ fðx1 ; x2 Þ : 0 6 x1 ; x2 6 2g,

8 1; > > > > > ðx2  x1 Þ; < 2  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v ð x1 ; x2 Þ ¼ > 1 2 2 þ 1 ; cos 4 p ð x 1  1:5Þ þ ðx2  0:5Þ > 2 > > > : 0;

x2  x1 P 0:5; 0 6 x2  x1 6 0:5; 1 ðx1  1:5Þ2 þ ðx2  0:5Þ2 6 16 ;

otherwise;

shown in Fig. 5. The derivative of this function is discontinuous on the straight lines x2  x1 ¼ 0:5 and x2  x1 ¼ 0 and on the circumferential line ðx1  1:5Þ2 þ ðx2  0:5Þ2 6 116. Besides, this function has local extrema and flat parts. This ‘‘rich’’ function is measured at 900 equidistant points in ½0; 2  ½0; 2. The resulting IIMLS surface is shown in Fig. 6 (a). It evidently shows the localized influence of the points of discontinuity, and the power of reconstructing the major

3124

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

Fig. 5. Graph of the non-smooth function.

Fig. 6. Graph of the IIMLS interpolated functions: (a) without truncation of influence domain and (b) with truncation of influence domain.

features of the approximated function. Besides, with clearly observations we can found that the IIMLS interpolated function is smooth on the lines x2  x1 ¼ 0:5 and x2  x1 ¼ 0. Fig. 7 (a) depicts the associated error v  v h . It can be found that the errors in the neighborhood of the lines x2  x1 ¼ 0:5 and x2  x1 ¼ 0 are considerable. This fact matches the smoothness of the approximated solution shown in Fig. 5. The reason for this phenomenon is that the above used IIMLS method is associated with interpolation of non-smooth functions. While the approximated function in this problem is continuous on the whole domain, its derivative is typically discontinuous across some lines. In the above computations, influence domains were allowed to extend across such discontinuous lines. Since the IIMLS shape functions and their derivatives generated with continuous basis functions are very smooth in the influence domain of an evaluation point, the interpolated function and its derivatives are continuous on these lines and thus some additional errors are inevitable.

Fig. 7. Error between true and interpolated surfaces: (a) without truncation of influence domain and (b) with truncation of influence domain.

3125

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

Fig. 8. Convergence for fitting of the non-smooth function.

To overcome this problem, the influence domain of an evaluation point can be restricted to the panel on which the approximated function is smooth. Namely, the influence domain of an evaluation point is truncated whenever it crosses a line or node where the derivative is discontinuous. This scheme is very beneficial to BIEs-based meshless methods, since the normal derivative of the primary unknown variable is typically discontinuous across a corner or an edge even if the primary variable is continuous. Such ideas have been successfully implemented in the BNM [12,13] and the GBNM n [14–17]. Now, the whole domain D is separated into four panels D1 ¼ fx2  x1 P 0:5g; D2 ¼ f0 6 x2  x1 6 0:5g, D3 ¼ ðx1  1:5Þ2 þ ðx2  0:5Þ2 6 116g and D4 ¼ D=ðD1 [ D2 [ D3 Þ. Then, the IIMLS method is performed on each panel independently. Fig. 6 (b) and Fig. 7 (b) give the corresponding IIMLS surface and error, respectively. Both figures clearly show that the errors in the neighborhood of the lines x2  x1 ¼ 0:5 and x2  x1 ¼ 0 disappear when the influence domain is truncated at discontinuous lines. The convergence of the usual scheme (without truncation) and the modified scheme (with truncation) is shown in Fig. 8, from which we can see that the modified scheme yields lower errors and has a higher convergence rate. 3. Interpolating boundary element-free method for three-dimensional potential problems 3.1. Boundary integral formulations Consider the following potential problem in X  R3 bounded by a piecewise smooth boundary C,

8 DuðxÞ ¼ 0; in X; > < ðxÞ; uðxÞ ¼ u on Cu ; > : @uðxÞ  q x ¼ ð Þ; on Cq ; qðxÞ :¼ @n ðxÞ

ð23Þ

 are the prescribed potential and normal flux, respectively, on the Dirichlet boundary Cu and on the Neumann  and q where u boundary Cq ¼ C  Cu , and n ¼ ðn1 ; n2 ; n3 Þ is the unit outward normal to C. Problem (23) can be used to formulate the general class of problems associated with the theory of gravitation, electrostatics, dielectrics, magnetostatics, fluid flow, flow in porous media and heat conduction [28]. For the potential u satisfying problem (23), from Green’s formula we have the following integral equation [35]:

uðxÞ ¼

Z

U ðx; yÞqðyÞdsy 

Z

C

T y ðx; yÞuðyÞdsy ;

x 2 X;

ð24Þ

C

where U ðx; yÞ ¼ 1=ð4pjx  yjÞ is the fundamental solution of the Laplace operator which is weak singular on the diagonal fðx; yÞ 2 C  C : x ¼ yg, and T y ðx; yÞ ¼ @U ðx; yÞ is the strongly singular fundamental solution. With the boundary conditions of problem (23), Eq. (24) can be further represented as

uðxÞ ¼

Z

U ðx; yÞqðyÞdsy þ

Cu

Z Cq

ðyÞdsy  U ðx; yÞq

Z

 ðyÞdsy  T y ðx; yÞu

Cu

Z

T y ðx; yÞuðyÞdsy ;

x 2 X:

ð25Þ

Cq

From Eq. (25) it is clear that if we want to recover the potential u in X we have firstly to find the yet unknown Cauchy data q on Cu and u on Cq . Sending in Eq. (24) x ! C, we have the following self-regular BIE [12]:

Z C

U ðx; yÞqðyÞdsy 

Z C

T y ðx; yÞ½uðyÞ  uðxÞdsy ¼ 0;

x 2 C:

ð26Þ

3126

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

Via first applying Eq. (26) on Cu and Cq independently, then inserting boundary conditions of problem (23), finally taking all unknown variables to the left side while all known variables to the right, one obtains:

Z

U ðx; yÞqðyÞdsy 

Z

Cu

T y ðx; yÞuðyÞdsy

ð27Þ

Cq

Z

¼

 ðy Þ  u  ðxÞdsy  T y ðx; yÞ½u

Z

Cu



 ðyÞ dsy ;  ðxÞ þ U ðx; yÞq T y ðx; yÞu

x 2 Cu ;

ð28Þ

Cq

and

Z

U ðx; yÞqðyÞdsy þ uðxÞ

Z

Cu

T y ðx; yÞdsy 

Cu

Z

¼

T y ðx; yÞ½uðyÞ  uðxÞdsy

ð29Þ

Cq

Z

 ðyÞdsy  T y ðx; yÞu

Z

Cu

ðyÞdsy ; U ðx; yÞq

x 2 Cq :

ð30Þ

Cq

3.2. The IIMLS scheme for the IBEFM Assume that the boundary face C is the union of piecewise smooth segments called panels. On each panel, curvilinear coordinates are defined as s ¼ ðs1 ; s2 Þ. For the actual calculations, s is taken to be the local relative coordinate of a boundary point with respect to a given evaluation point. Otherwise, a loss of accuracy follows from the absolute values of s being too large with respect to one. Therefore, the value of s for an evaluation point is always 0. In order to avoid the discontinuity at edges and corners, the influence domain of a boundary point is restricted to the panel to which it belongs. In this scheme, collocation nodes are located inside panels, and no node is placed on the edge or at the corner of the boundary. Let fxi gNi¼1 denote an arbitrarily chosen set of N boundary nodes xi 2 C. To simplify the representation, we further assume u that the first N u boundary nodes fxi gNi¼1 belong to Cu and the rest N  N u boundary nodes fxi gNi¼Nu þ1 belong to Cq . Based on fxi gNi¼1 alone, the IIMLS scheme for the potential u and its normal gradient q ¼ @u=@n can be defined as

qðxðsÞÞ 

Nu X

Ui ðsÞqi ;

xðsÞ 2 Cu ;

ð31Þ

i¼1

uðxðsÞÞ 

N X

Ui ðsÞui ;

xðsÞ 2 Cq ;

ð32Þ

i¼N u þ1

where ui and qi are, respectively, the actual nodal values uðxi Þ and qðxi Þ. Ui ðsÞ is the contribution from the node xi to the evaluation point xðsÞ. In the influence domain of the ith node, Ui ðsÞ is the IIMLS shape function given by Eq. (17) and is not equal to zero. 3.3. Discretization of BIEs u Inserting Eqs. (31) and (32) into Eq. (27) for boundary nodes fxk gNk¼1  Cu leads to

Z Nu X qi Cu

i¼1

N X

U ðxk ; yÞUi ðyÞdsy 

uj

j¼N u þ1

Z

T y ðxk ; yÞUj ðyÞdsy ¼

Cq

Z

 ðy Þ  u  ðxk Þdsy T y ðxk ; yÞ½u

Cu



Z

  ðyÞ dsy :  ðxk Þ þ U ðxk ; yÞq T y ðxk ; yÞu

ð33Þ

Cq

While inserting Eqs. (31) and (32) into Eq. (29) for boundary nodes fxm gNm¼Nu þ1  Cq yields

Z Nu X qi i¼1

¼

Z

N X

U ðxm ; yÞUi ðyÞdsy þ

Cu

 ðyÞdsy  T y ðx; yÞu

Cu

Z

(

uj Uj ðxm Þ

j¼N u þ1

Z

T y ðxm ; yÞdsy  Cu

Z

  T y ðxm ; yÞ Uj ðyÞ  Uj ðxm Þ dsy

)

Cq

ðyÞdsy : U ðx; yÞq

ð34Þ

Cq

Eqs. (33) and (34) result in the following N  N linear system of equations



A11

A12

A21

A22

 q u

 ¼

fu fq

;

ð35Þ

3127

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

 T where q ¼ q1 ; q2 ; . . . ; qNu is the vector of N u unknown nodal values qi ; u ¼ ½uNu þ1 ; uNu þ2 ; . . . ; uN T is the vector of N  N u unknown nodal values uj . The block matrices used in Eq. (35) are given by

A11 ½k; i ¼

Z

U ðxk ; yÞUi ðyÞdsy ;

ð36Þ

Cu

Z

A12 ½k; j ¼ 

T y ðxk ; yÞUj ðyÞdsy ;

ð37Þ

Cq

Z

A21 ½m; i ¼

U ðxm ; yÞUi ðyÞdsy ;

ð38Þ

Cu

A22 ½m; j ¼ Uj ðxm Þ

Z

T y ðxm ; yÞdsy 

Cu

Z

  T y ðxm ; yÞ Uj ðyÞ  Uj ðxm Þ dsy ;

ð39Þ

Cq

for all k; i ¼ 1; . . . ; N u and m; j ¼ N u þ 1; . . . ; N. The components of the right-hand side are given by

f u ½k ¼

Z

 ðy Þ  u  ðxk Þdsy  T y ðxk ; yÞ½u

Cu

f q ½m ¼

Z Cu

Z



 ðyÞ dsy ;  ðxk Þ þ U ðxk ; yÞq T y ðxk ; yÞu

ð40Þ

Cq

 ðyÞdsy  T y ðx; yÞu

Z

ðyÞdsy ; U ðx; yÞq

ð41Þ

Cq

for all k ¼ 1; . . . ; N u and m ¼ N u þ 1; . . . ; N. Eqs. (36)–(41) have integrations over the boundary. These integrations can be calculated numerically by employing a cell structure as in the EFG and the BNM. The integrands in Eqs. (37), (38) and (40) are always regular and can be evaluated by usual Gaussian quadrature. The integrands in Eqs. (36), (39) and (41) consist of regular and weakly singular functions. The regular integral can be computed in the same manner as above, while the weakly singular integral may be handled using regularization procedures [12,13]. Once the unknown vectors q and u are found from Eq. (35), the yet unknown Cauchy data q on Cu and u on Cq could be computed using Eqs. (31) and (32), respectively. Then, the potential and its derivatives at an internal point x 2 X may be computed using Eq. (25). 3.4. Comparison between the IBEFM, the BNM and the BEM A comparison between the IBEFM, the BNM and the BEM is summarized concisely in Table 1. All methods are based on BIEs. The difference is the way in which the shape function is formulated. The IBEFM is a boundary type meshless method, whereas the BEM is a boundary type numerical method based on boundary elements. As other meshless methods (e.g. the EFG, the MLPG and the BNM), the interpolation procedure in the IBEFM is based only on a cluster of boundary nodes as discussed above. But in the BEM, interpolants of boundary variables are related to the geometry of the elements. Compared with the BEM, the IBEFM does not require an element structure and thus simplifies mesh generation and problems which require modification of the mesh, such as adaptive refinement. Both the IBEFM and the BNM are boundary type meshless methods. The shape function of the BNM generated using the MLS scheme lacks the delta function property, thus extra effort is needed in the BNM to enforce boundary conditions. This extra effort is avoided in the IBEFM. Besides, the BNM is an indirect numerical method in which the basic unknown quantity is approximations of nodal variables, while the IBEFM is a direct numerical method in which the basic unknown quantity is real solutions of nodal variables. From Table 1, we can find that the number of both system equations and unknowns in the BNM is double of those in the IBEFM. This is because N additional equations have to be employed in the BNM to impose boundary conditions. As in the BEM, system matrix of the BNM and the IBEFM is dense and unsymmetrical, requiring

Table 1 Comparison between the IBEFM, the BNM and the BEM.

Mesh Interpolation Interpolation domain based on Delta property of shape function Implementation of boundary conditions Number of system equations for N boundary nodes Unknowns of system equations Number of unknowns of system equations

IBEFM

BNM

BEM

No IIMLS Distributed boundary nodes Yes Easy N Actual nodal values N

No MLS Distributed boundary nodes No Difficult 2N Approximate nodal values 2N

Yes Polynomial Pre-defined elements Yes Easy N Actual nodal values N

3128

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

Fig. 9. Errors and convergence rates for 3D Dirichlet problem on a cube: (a) u, (b) @u=@x1 , (c) @u=@x2 and (d) @u=@x3 .





O N 2 memory and O N 3 operations. Therefore, the memory required for storing the system equations in the IBEFM is only a quarter of that in the BNM, while the CPU time required for solving the system equations in the IBEFM is only eighth that in the BNM. 4. Numerical experiments 4.1. Dirichlet and mixed problems on a cube The first example embraces a variety of Dirichlet and mixed problems on a cube. The cube faces are x1 ¼ 1; x2 ¼ 1 and x3 ¼ 1, respectively. The following solutions, which are taken from Refs. [12,13], are employed as analytical solutions: (i) Linear solution:

u ¼ x1 þ x2 þ x3 ;

ð42Þ

(ii) Quadratic solution-1:

u ¼ x1 x2 þ x2 x3 þ x3 x1 ;

ð43Þ

(iii) Quadratic solution-2:

u ¼ 2x21 þ x22 þ x23 ;

ð44Þ

(iv) Cubic solution:

u ¼ x31 þ x32 þ x33  3x21 x2  3x22 x3  3x23 x1 :

ð45Þ

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

3129

Fig. 10. Results of u along the x1 -axis for 3D mixed problem on a cube.

First, Dirichlet problems are solved for which the essential boundary conditions are imposed on all faces corresponding to the analytical solutions. Fig. 9 plots the log–log plot of errors of u and its derivatives @u=@x1 with respect to the number of boundary nodes. The results are obtained using five sets of nodes: (a) 48 nodes, (b) 96 nodes, (c) 192 nodes, (d) 384 nodes and (e) 768 nodes. It can be found that the experimental convergence rate is very high and approximate to 2.5, and the rates are almost identical for all analytical solutions. Besides, u and its derivatives have almost the same experimental convergence rates. Next, boundary value problems with mixed boundary conditions are also considered. Dirichlet boundary conditions are imposed on faces x3 ¼ 1 and Neumann boundary conditions on all other faces. The analytical and numerical solutions for the potential u and its derivative @u=@x1 are exhibited in Figs. 10 and 11. The linear, quadratic-2 and cubic solutions are used here. The results are obtained using 32 nodes on each face. From these figures, it can be observed that the numerical solutions are in excellent agreement with the analytical solutions. As in Refs. [12,13], the IBEFM has also been applied on a more challenging problem which cannot be described via polynomial approximations. The analytical solution of this problem is

u ¼ sinh

p x1 2

px2 px3 sin pffiffiffi sin pffiffiffi : 2 2 2 2

Dirichlet boundary conditions are imposed on the surface of the cube. Fig. 12 depicts the log–log plot of errors with respect to the number of boundary nodes. Figs. 13 and 14 exhibit CPU times required for constructing and solving system equations, respectively. In these figures, the results of the BNM are also given for comparison. In this analysis, we employed three sets of nodes: (a) 192 nodes, (b) 768 nodes and (c) 3072 nodes. Fig. 12 shows that the solution accuracy of the IBEFM is higher than that of the BNM, while Figs. 13 and 14 show that the IBEFM requires less CPU times than the BNM for both constructing and

solving system equations. Moreover, the two methods nearly require O N 2 operations for constructing coefficient matrices

and O N 3 operations for solving system equations. Consequently, the CPU time for solving the system equations in the BNM is about eightfold that in the IBEFM. 4.2. Problems on a horn domain The current IBEFM has also been tried on more complicated geometries. In this example, the boundary surface is a horn. Fig. 15 depicts the sketched geometric configuration and a discretization consisting of 288 integration cells with one node

Fig. 11. Results of @u=@x1 along the x1 -axis for 3D mixed problem on a cube.

3130

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

Fig. 12. Comparison of errors for a Dirichlet problem on a cube: (a) u and (b) @u=@x1 .

Fig. 13. Comparison of CPU times for constructing system equations for a problem on a cube.

Fig. 14. Comparison of CPU times for solving system equations for a problem on a cube.

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

3131

Fig. 15. Schematic diagram and integration cells for a horn domain (the dot represents boundary nodes).

per cell. To assess the accuracy of the IBEFM for this problem, the following oscillatory harmonic function is considered as the analytical solution,

uðx1 ; x2 ; x3 Þ ¼ ex1 sin x3 þ ex3 cos x2 : Dirichlet boundary conditions corresponding to the analytical solution are prescribed on the side surface of the horn and Neumann boundary conditions on the top and bottom surfaces. Fig. 16 gives a comparison between the present numerical results with the analytical solutions for u and its derivatives along the arc given by the formulas x1 ¼ 9ð1 þ sin hÞ; x2 ¼ 0; x3 ¼ 9ð1 þ cos hÞ; h 2 ½p; 0:5p. This figure shows that the numerical solutions are seen to capture the behavior of the analytical solutions very well.

Fig. 16. Results of (a) u, (b) @u=@x1 , (c) @u=@x2 and (d) @u=@x3 for the problem on a horn domain.

3132

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

Fig. 17. Schematic diagram and integration cells for a cylindrical tube (the dot represents boundary nodes).

4.3. Problems on a cylindrical tube Finally we consider a mixed boundary value problem for the Laplace’s equation in a cylindrical tube X centered at the origin Oð0; 0; 0Þ such that



X ¼ ðx1 ; x2 ; x3 Þ 2 R3 : r <

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x21 þ x22 < R; H < x3 < H :

Fig. 17 depicts the sketched geometric configuration and a discretization consisting of 864 integration cells with one node per cell. The cubic polynomial, Eq. (26), is employed as the analytical solution. Dirichlet conditions are imposed on the top qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n o face r 6 x21 þ x22 6 R; x3 ¼ H and Neumann conditions on the remaining surface components of the whole boundary surface. The parameters are taken as r ¼ 1, R ¼ 2 and H ¼ 4 in the computation.

Fig. 18. Results of (a) u, (b) @u=@x1 , (c) @u=@x2 and (d) @u=@x3 for the problem on a cylindrical tube.

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

3133

Fig. 18 exhibits a comparison between numerical results with analytical solutions for u and its derivatives of some boundary points. Those points are located along the line segment from ð2; 0; 4Þ to ð0; 2; 4Þ. Again, we can observe that the numerical solutions are in excellent agreement with the analytical solutions. 5. Conclusions Based on combining BIEs with the IIMLS method, an interpolating boundary element-free method (IBEFM) is presented in this paper. The IBEFM takes the advantages of both BIEs in dimension reduction and the IIMLS method in elements elimination. The properties of the IIMLS shape function are discussed in detail. Curve and surface fitting capability of the IIMLS method is evaluated through one- and two-dimensional functions. Convergence studies indicate that the IIMLS method possesses a remarkable accuracy and has high rates of convergence. These convergence rates are higher than those from the MLS approximation. Besides, numerical examples are presented for three-dimensional potential problems. The numerical results have verified the accuracy, convergence and effectiveness of the IBEFM. The errors and CPU times were compared to those obtained by the BNM. These comparisons show that the IBEFM provides monotonic convergence and high accuracy results with low computational expenses. The present IBEFM is extensible to problems in real world systems such as solid mechanics, fluid mechanics, wave propagation, electromagnetic problems, and so on. Acknowledgements This work was supported by National Natural Science Foundation of China under Grant Nos. 11471063 and 11101454, the Natural Science Foundation Project of CQ CSTC under Grant Nos. cstc2014jcyjA00005 and cstc2013jcyjA30001, and the Program of Chongqing Innovation Team Project in University under Grant No. KJTD201308. The author is very grateful to the anonymous reviewers for their invaluable comments and suggestions which have greatly improved the quality and the presentation of this paper. References [1] P. Lancaster, K. Salkauskas, Surface generated by moving least squares methods, Math. Comput. 37 (1981) 141–158. [2] X. Li, J. Zhu, A Galerkin boundary node method and its convergence analysis, J. Comput. Appl. Math. 230 (2009) 314–328. [3] X. Li, Meshless Galerkin algorithms for boundary integral equations with moving least square approximations, Appl. Numer. Math. 61 (2011) 1237– 1256. [4] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl, Meshless methods: an overview and recent developments, Comput. Methods Appl. Mech. Eng. 139 (1996) 3–47. [5] C.A. Duarte, J.T. Oden, H-p clouds—An h-p meshless method, Numer. Methods Partial Diff. Eq. 12 (1996) 673–705. [6] S.N. Atluri, The Meshless Method (MLPG) for Domain & BIE Discretizations, Tech. Science Press, California, 2004. [7] S. Li, W.K. Liu, Moving least square reproducing kernel method (II) Fourier analysis, Comput. Methods Appl. Mech. Eng. 139 (1996) 159–193. [8] W.K. Liu, S. Li, T. Belytschko, Moving least square reproducing kernel method (I) methodology and convergence, Comput. Methods Appl. Mech. Eng. 143 (1997) 113–154. [9] S. Li, W.K. Liu, Reproducing kernel hierarchical partition of unity part I: formulations, Int. J. Numer. Methods Eng. 45 (1999) 251–288. [10] S. Li, W.K. Liu, Reproducing kernel hierarchical partition of unity part II: applications, Int. J. Numer. Methods Eng. 45 (1999) 289–317. [11] Y.X. Mukherjee, S. Mukherjee, The boundary node method for potential problems, Int. J. Numer. Methods Eng. 40 (1997) 797–815. [12] S. Mukherjee, Y.X. Mukherjee, Boundary Methods: Elements, Contours, and Nodes, CRC Press, Boca Raton, 2005. [13] M.K. Chati, S. Mukherjee, The boundary node method for three-dimensional problems in potential theory, Int. J. Numer. Methods Eng. 47 (2000) 1523– 1547. [14] X. Li, J. Zhu, A Galerkin boundary node method for biharmonic problems, Eng. Anal. Boundary Elem. 33 (2009) 858–865. [15] X. Li, Meshless analysis of two-dimensional Stokes flows with the Galerkin boundary node method, Eng. Anal. Boundary Elem. 34 (2010) 79–91. [16] X. Li, Adaptive meshless Galerkin boundary node methods for hypersingular integral equations, Appl. Math. Model. 36 (2012) 4952–4970. [17] X. Li, The meshless Galerkin boundary node method for Stokes problems in three dimensions, Int. J. Numer. Methods Eng. 88 (2011) 442–472. [18] K.M. Liew, Y. Cheng, S. Kitipornchai, Boundary element-free method (BEFM) and its application to two-dimensional elasticity problems, Int. J. Numer. Methods Eng. 65 (2006) 1310–1332. [19] M. Peng, Y. Cheng, A boundary element-free method (BEFM) for two-dimensional potential problems, Eng. Anal. Boundary Elem. 33 (2009) 77–82. [20] H. Ren, Y. Cheng, W. Zhang, An improved boundary element-free method (IBEFM) for two-dimensional potential problems, Chin. Phys. B 18 (2009) 4065–4073. [21] H. Ren, Y. Cheng, W. Zhang, An interpolating boundary element-free method (IBEFM) for elasticity problems, Sci. Chin. Ser. G Phys. Mech. Astron. 53 (2010) 758–766. [22] I. Kaljevic, S. Saigal, An improved element free Galerkin formulation, Int. J. Numer. Methods Eng. 40 (1997) 2953–2974. [23] H. Ren, Y. Cheng, The interpolating element-free Galerkin (IEFG) method for two-dimensional elasticity problems, Int. J. Appl. Mech. 3 (2011) 735–758. [24] H. Ren, Y. Cheng, The interpolating element-free Galerkin (IEFG) method for two-dimensional potential problems, Eng. Anal. Boundary Elem. 36 (2012) 873–880. [25] H. Netuzhylov, Enforcement of boundary conditions in meshfree methods using interpolating moving least squares, Eng. Anal. Boundary Elem. 32 (2008) 512–516. [26] J. Wang, F. Sun, Y. Cheng, An improved interpolating element-free Galerkin method with nonsingular weight function for two-dimensional potential problems, Chin. Phys. B 21 (9) (2012) 090204. [27] J. Wang, J. Wang, F. Sun, Y. Cheng, An interpolating boundary element-free method with nonsingular weight function for two-dimensional potential problems, Int. J. Comput. Methods 10 (2013) 1350043. [28] A.P.S. Selvadurai, Partial Differential Equations in Mechanics, vol.1, Springer, Berlin, 2000. [29] G.R. Liu, Y.T. Gu, A point interpolation method for two-dimensional solids, Int. J. Numer. Methods Eng. 50 (2001) 937–951. [30] X. Li, A meshless method based on boundary integral equations and radial basis functions for biharmonic-type problems, Appl. Math. Model. 35 (2011) 737–751. [31] Y.T. Gu, G.R. Liu, A boundary point interpolation method for stress analysis of solids, Comput. Mech. 28 (2002) 47–54.

3134

X. Li / Applied Mathematical Modelling 39 (2015) 3116–3134

[32] X. Li, J. Zhu, S. Zhang, A hybrid radial boundary node method based on radial basis point interpolation, Eng. Anal. Boundary Elem. 33 (2009) 1273– 1283. [33] X. Li, Numerical solution of solid mechanics problems using a boundary-only and truly meshless method, Math. Prob. Eng. 2012 (2012). Article ID: 298903. [34] R. Scitovski, Š. Ungar, D. Kukic, Approximating surface by moving total least squares method, Appl. Math. Comput. 93 (1998) 219–232. [35] J. Zhu, Z. Yuan, Boundary Element Analysis, Science Press, Beijing, 2009.