h-Adaptive radial basis function finite difference method for linear elasticity problems

In this research work, the radial basis function finite difference method (RBF-FD) is further developed to solve one- and two-dimensional boundary value problems in linear elasticity. The related differentiation weights are generated by using the extended version of the RBF utilizing a polynomial basis. The type of the RBF is restricted to polyharmonic splines (PHS), i.e., a combination of the odd m-order PHS ϕ(r)=rm\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi (r)=r^m$$\end{document} with additional polynomials up to degree p will serve as the basis. Furthermore, a new residual-based adaptive point-cloud refinement algorithm will be presented and its numerical performance will be demonstrated. The computational efficiency of the PHS RBF-FD approach is tested by means of the relative errors measured in ℓ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _2$$\end{document}-norm on several representative benchmark problems with smooth and non-smooth solutions, using h-adaptive, uniform, and quasi-uniform point-cloud refinement.


Introduction
In the last two decades, the application of the RBF technique as a mesh-free scheme has become popular for the numerical solution of partial differential equations (PDE). The fundamental concept of RBF-based numerical approximation of differential operators originates from the research works [41,42], and a large number of papers have since contributed to the continuous development of the RBF-FD method [8,17,19,34,36,37,43] and dealt with the approximative solution of numerous BVPs and initial-BVPs appearing in science and engineering [6,25,27,38]. The main advantage of the advanced FD methods based on the RBFs is that they are well-suited for arbitrarily-scattered (non-structured) point layouts. Furthermore the stencil topology can be chosen in different ways, thereby avoiding the fixed rectangular  [21].
During the application of the RBF-FD scheme, the socalled differential weights are locally computed at a certain evaluation point within a stencil, using an RBF-based approximation. In this case, the stencil is constructed by a fixed number of points that are closest to the evaluation point under consideration. Accordingly, we start off by computing the differential weights at each point of the point-cloud generated on the considered computational domain, before storing them for the global assembling procedure-which is the next important step in solving the PDEs. Thus, the weight computations are performed as the pre-processing step of the solution algorithm. Since this computation can be carried out independently for all evaluation points, this part of the RBF-FD method is well parallelizable.
The RBFs can be divided into two main groups: piecewiseand infinitely smooth RBFs. The latter ones-such as the Gaussian-, the multiquadratic and inverse quadratic and inverse multiquadratic RBFs-provide spectral convergence behavior for smooth BVPs. However, they are not independent of a so-called shape parameter, which controls the gradient of the RBFs, thereby having significant effect on the condition number of the RBF matrices and thus on the solution accuracy. The optimal value of the shape parameter has to be found for all evaluation points during the weight computation in order to achieve numerical stability and accuracy [5,7,10,12,22]. The RBF coefficient matrix conditioning problem can be overcome by a stabilization algorithm using small shape parameter values [20,26] or by applying RBFs with a hybrid kernel [29].
The piecewise smooth RBFs (such as the PHS) have a great advantage over the infinitely smooth RBFs, since they are not dependent on a shape parameter and since their usage can yield higher-order accuracy-though only with algebraic h-convergence characteristics. Additionally, the piecewise smooth RBFs perform well for the approximation of the differential operators near to the boundaries where the stencil geometry becomes half-sided, i.e., the deterioration of the numerical approximation accuracy due to the Runge phenomenon can be avoided, see [1][2][3][4].
In accordance with the above-mentioned beneficial properties of the piecewise smooth RBFs, in this paper, the PHS RBFs supplemented with polynomial functions will be introduced as the basis function space in order to develop an advanced FD scheme for the numerical solution of 1D and 2D BVPs in linear elasticity. The robust convergence behavior of the FD method based on this type of augmented basis function spaces has been verified in [1,2,4,14]. Furthermore, it has been experienced in numerical studies presented in [4,14] that the highest polynomial degree specifies the convergence rate, while the PHS RBF order has a significant influence only on the achievable accuracy. One order in the convergence rate is lost in the case of the numerical approximation of the first (spatial) derivatives, while two orders are lost for approximating the second derivative. Thus, the convergence rate depends on the order of the differential operator [4,5,13,14,30].
Interestingly, the increase in the stencil size does not lead to better accuracy and faster convergence when a fixed higher-degree polynomial basis is used. However, if lowerdegree polynomials are applied and their degree is kept fixed, increasing the stencil size can yield improved accuracy while the convergence order remains unchanged [2]. This type of independence on the stencil size motivated the construction of the p-adaptive PHS RBF-FD method. Here the desired accuracy is achieved by locally increasing the polynomial degree and adapting the required stencil size, see [30], for which the non-uniform point-cloud is generated by the code NodeLab [28].
Adaptive strategies such as h-, p-and hp-refinement are popular for finite element methods (FEMs), see [11,35] for example. The basic idea of the adaptive schemes is to introduce an error indicator function that guides the non-uniform h-and/or p-refinement until the prescribed stopping criteria are met. Although the FEM-based adaptive strategies are more prevalent than the FD method-based ones, the construction of the adaptive point refinement procedure is much easier for RBF-FD schemes than the implementation of the adaptive mesh refinement process in the case of FEMs. Residual errorand Zienkiewicz-Zhu type indicator-based adaptive RBF-FD algorithms are also developed using the Gaussian RBF, for example in [23,24,31]. It is worth emphasizing here that, in the case of h-adaptive point refinement, the DistMesh code is one of the most suitable for fast generation of 2D point distributions with variable density [15,31,33].
In accordance with the above findings, the aim of this paper is to develop a residual-based h-adaptive point refinement algorithm for 1D BVPs, which is based on the odd m-order PHS supplemented with additional polynomials. Firstly, the computational efficiency of the approach will be investigated on smooth and non-smooth 1D BVPs, using the uniform-and the adaptive point-cloud refinement schemes. Secondly, the h-convergence behavior of the PHS RBF-FD method complemented with polynomials will be tested on smooth and non-smooth 2D BVPs in linear elasticity, applying uniform, quasi-uniform and h-adaptive point distributions.
The paper is organized as follows. In Sect. 2, the basic system of second-order differential equations and the related Dirichlet-and Neumann-type boundary conditions are summarized for 2D linear elasticity problems. Section 3 briefly presents the procedure to evaluate the local RBF matrix and to compute the related differentiation weights from the local linear algebraic system of equations using PHS combined with polynomials. In Sect. 4, a comprehensive overview of the applied stencil geometry is given, discussing how to treat the Neumann-type boundary conditions with the aid of ghost points. Furthermore, a new adaptive point-cloud refinement strategy is presented. In Sect. 5, the h-convergence behavior of the PHS RBF-FD method is studied for 1D BVPs, applying uniform as well as the proposed adaptive pointcloud refinement. Afterwards, the PHS RBF-FD method is also investigated for 2D BVPs of linear elasticity-applying the uniform, the quasi-uniform, and the newly-developed h-adaptive point-cloud refinements. Section 6 summarizes the results and suggests possible extensions of the presented approach.

Basic system of differential equations in plane elasticity
The second-order partial differential equations for the 2D elastostatic problem are expressed in a Cartesian coordinate frame as where u x , u y , and f x , f y are the coordinates of the displacement-and body force density vectors u = u x i +u y j and f = f x i + f y j . Further, μ and λ denote the Lamé coefficients, while L as well as L x x , L yy , and L xy are, respectively, the Laplacian as well as the pure and mixed 2D second-order partial differential operators-defined as The coefficient matrix in Eq. (1) can be rearranged in the form The Neumann-type boundary conditions prescribed on the surface part S p ⊂ ∂ are expressed in terms of the firstorder partial derivatives of the displacements u x and u y : using the operator-notations Here, ∂ = S p ∪ S u is the boundary of the body for which S p ∩ S u = ∅ is valid, while n x and n y denote the coordinates of the outward unit normal vector n = n x i + n y j to ∂ . The Dirichlet-type boundary condition u =ũ =ũ x i +ũ y j on S u is prescribed on the remaining boundary part of the solid body S u ⊂ ∂ . The Lamé coefficient λ corresponds to if plane stress conditions are considered. Here, λ converges to ∞, as ν tends to 1 2 in the plain strain case, i.e., this model problem will be prone to volumetric locking effects.

Computing the differentiation weights using the augmented PHS
In this section, our aim is to approximate the first-and secondorder partial derivative operators L introduced in Eqs. (2) and (5) at the evaluation point r e . The derivative Lu at r e is approximated by a linear combination of the function values u i = u(r i ) considered at the distinct points r i ∈ R d along with the evaluation point r e , where i = 1, . . . , n and d is the space dimension. These points are used to construct the local stencil. Accordingly, we seek the weights w i as the coefficients of u i in such a way that in which n and w i are called the stencil size and the differentiation weights, respectively. In order to determine w i , we start with a linear combination along with the constraints where P j (r i ) denotes multivariate polynomials up to degree p and φ( r − r i ) are the RBFs localized at the points r i , [1,13,16,30]. Here, r − r i = r − r i 2 denotes the Euclidean norm on R d . For the mathematical background of the related quadratic constrained minimization problem, we refer to [2,4].
To derive the equation for the differentiation weights we apply the interpolant (7) to the function values u i = u(r i ) and consider the p+d p constraints (8). Next, we use again the interpolant (7) but now for the approximation of the derivative [Lu(r)] r=r e at the evaluation point r e . Then we eliminate the function values u i with the assumption that k i and γ i are non-zero coefficients arriving at the following linear equation system which can be solved for the weights w i and v j , for details see [1,14]. In Eq. (9) represent the symmetric RBF matrix K ∈ R n×n and the are the differentiation vectors ∈ R n and q ∈ R ( p+d p ) , as well as . Please note that the weights v are omitted, for details see [1,14,16].
The widely-used RBFs are classified into two main groups: the infinitely and piecewise smooth RBFs, see the tables in [1,13,16]. The infinitely smooth RBFs depend on the shape parameter , influencing the gradient of the RBFs and thereby controlling the condition number of the coefficient matrix in Eq. (9) and thus the accuracy of the result for the differentiation weights. On the contrary, the piecewise smooth PHS is independent of and provides higher convergence rates under node refinement, without the use of a stable algorithm to invert the RBF matrix, see [1][2][3][4]18]. This is the main reason and motivation why the PHS is chosen for the numerical computation of the weights, supplemented with a polynomial basis for the sake of further stabilization and faster convergence.
The RBF-FD methods are well-applicable if the coefficient matrix appearing in Eq. (9) exhibits non-singular behavior, i.e., the system (9) has a unique solution for the choice of the odd m-order PHS in which r e = x e i + y e j represents the Cartesian coordinates of the evaluation point of the stencil, supplemented with polynomials up to degree p where If, for instance, m = 3 is chosen, the degree of the polynomial basis has to be at least 1-or degree p has to be at least 2 for m = 5, see [1]. The stencil size n is chosen to be approximately twice the number of the polynomial basis functions p+d p , which is p + 1 in 1D and ( p + 2)( p + 1)/2 in the 2D case. The increase in stencil size n has only small or no influence on the accuracy and convergence rate for a fixed polynomial degree p, see [4].

Point cloud and stencil geometry
In what follows, the applied point-cloud will be structured as Cartesian or, in other words, as a rectangular point distribution for the convergence studies of the uniform h-refinement, see Figs. 1 and 7. These Cartesian arrangements are similar to the standard grid-like FD schemes, but the RBF interpolant defined in Eqs. (7)-(8) is of course quite different. Next to the boundaries, the stencil becomes half-sided and non-symmetric. Although the evaluation point comes closer to the boundaries in this case the stencils retain the rectangular shape. However, inside the domain, the stencils remain symmetric and their evaluation point is localized at their center, see Figs. 1 and 7.
Next, we consider Neumann-type boundary conditions (BC). Let us assume that the Neumann-type BC (4) and its 1D version (15), respectively, are prescribed on the right boundary x = 2 and on the top boundary y = 1 for the 2D case ( Fig. 7), as well as at the right boundary x = L for the 1D case ( Fig. 1).
In order to satisfy these BCs at a given boundary point, a new point has to be added to the point-cloud outside the domain as a neighbor of the boundary point under consideration. These external points are called ghost points. They allow for the addition of constraint equations to the global system of equations arising from the numerical approximation of the 2D PDE (1) or its 1D version (13). The subsidiary equations are required for the enforcement of the Neumann-type BC (4) (or (15) for the 1D case) prescribed at the corresponding boundary points, see [14].
As a result of the latter procedure, new unknown function values have to be introduced at the external points for u in 1D and for u x and u y in 2D. After the global assembly process, the resulting linear algebraic system of equations is simultaneously solved for all of the unknown function values located at the external, internal, and boundary points. The total number of the points N in the point-cloud is given as the sum of the number of the external, internal, and boundary points, N E , N I , and N B : As described above, the number of the extra constraint equations associated with the enforcement of the Neumann-type BC is equal to d · N E , where d denotes the space dimension, see [14].

h-adaptive point-cloud refinement strategy
In this subsection, we propose a new adaptive point-cloud refinement scheme for the PHS-RBF method to solve 1D and 2D problems of linear elasticity. In what follows, the h-adaptive algorithm will be demonstrated for 1D problems. This type of BVP can be considered as the 1D version of the problem presented in Sect. 2. The basic equation and the boundary conditions can be stated as follows, see [11,39]: Here, f (x) represents as a distributed line load, A the crosssectional area, E the Young's modulus, and L the length of the bar-whileũ andF are the prescribed values of the axial displacement at x = 0 and the axial force at x = L, respectively. The point refinement procedure is guided by an error indicator based on the residual of the second-order differential equation (13), see Algorithm 1. At the very beginning of the numerical process, we set the initial number of points N 0 , the number of the iteration steps, the order of the PHS m, and the degree of the supplementary polynomials p, as well as the material parameters and the loading values. Next, the structured point-cloud x 0 is generated.
In each iteration step, the problem has to be solved. To this end, the weights for the approximation of L x and L x x at the internal and boundary points are computed by solving the local system (9) and assembling the global first-and second-order differential operator matrices D x , D x x and the global equation system AE D x x u = f. Next, Dirichlet-type boundary conditions (14) are applied in strong form, and Neumann-type boundary conditions (15) as new constraint equation are considered. Then, the resulting system is solved in order to determine u ∈ R N . Afterwards, the column matrix of the axial force F is computed as AE D x u. By means of the first-order differential-operator matrix D x , we compute the numerical approximation of the right-hand side of Eq. (13) asf = D x F. In this way, the residual is expressed in the form In a next step, the absolute value of the residual is checked for two neighboring points within a 'for' loop over all the internal and boundary points. A new point is inserted between the two neighboring points in question if the refinement criterion (17) holds true for one of the two investigated points. Here, η max is the maximum norm of the column matrix of the global residual as whereas N is the total number of points and β ∈ [0, 1] denotes the threshold parameter. The choice of this tolerance value controls how many points are added to the point-cloud in each iteration step. The h-adaptive refinement terminates if it meets the stopping criterion, which can be set for both the number of the iterations or the maximum of the residual. During our numerical tests, the value β will be set to 0.01 for 1D BVPs and to 0.1 for 2D BVP. Algorithm 1 summarizes the h-adaptive point-cloud refinement scheme.

Numerical experiments
In this section, the performance of the PHS-based RBF-FD method augmented with a polynomial basis will be analyzed for both smooth and non-smooth BVPs of 1D and 2D linear elasticity. During the computational tests, the convergence of the relative error measured in 2 -norm of the displacements and stresses will be studied. The number of degrees of freedom (denoted as DOF in the figure captions) is defined as the total number of points N of the point-cloud, ignoring the constraints related to the Dirichlet BC. Additionally, the influence of the polynomial degree p on the rate of h-convergence will be investigated in each problem. In order to avoid the singular behavior of PHS when x = x e and y = y e for a second-order operator, the PHS order m has to be chosen high enough. In what follows we set m to 3 for 1D and m = 5 for 2D BVPs. Therefore, the smallest possible polynomial degree p, on account of inequality (11), is 1 in 1D and 2 in 2D.

One-dimensional boundary value problems
In this subsection, the computational performance of the PHS-based RBF-FD method will be analyzed for the 1D problem presented above. For the sake of simplicity, the geometric parameters A and L, as well as Young's modulus E, will be assumed to be constant. In view of this, the values AE = 1 N and L = 1 m will be assumed in all of the considered 1D benchmark problems.
During the h-convergence study, the domain will be refined equidistantly in 10 steps for smooth and non-smooth Compute n = 2 p + 5 stencil size 3: Generate x 0 , using N 0 initial uniform point-cloud 4: Set iterstep maximum number of iterations 5: Set β tolerance value for the error 6: Set iter = 1 initial value of the iteration number 7: repeat 8: Compute the weights for L x and L xx at the internal and boundary points 9: AE D xx u = f assemble the global equation system Eq. (13) 10: Apply the BCs ũ andF according to Eqs. (14)-(15) 11: Solve the resulting linear algebraic equation system for the column matrix u 12: F = AE D x u compute the column matrix of the internal force F 13:f = D x F compute the right-hand side column matrixf 14: η =f − f compute the column matrix of the residual 15: η max = max |η i | determine the maximum of the residual according to Eq. (18) 16: N new = 0 set the initial number of the newly-inserted points to zero 17: for i = 1 to N (loop over all the internal and boundary points) do 18: Check whether a new point needs to be put between two neighboring points 19: iter = iter + 1 increment the iteration number 28: until ( (η max /N < β) or (iter > iterstep) ) 29: end procedure (such as singular-and shock) BVPs, while the degree of the polynomial basis p is fixed in each step. In order to investigate the effect of the polynomial degree p on the rate of the h-convergence, p will be varied from 1 to 7 and 2 to 5, respectively, for smooth and non-smooth problems. Therefore, the odd stencil size n has to be adjusted to the actual p-value as n = 2 p + 5.
The efficiency of the newly developed adaptive point refinement strategy will also be investigated by comparing its convergence behavior to the numerical results obtained using equidistant point sets.
The symmetric stencil geometry applied within the domain and the half-sided, non-symmetric stencil shapes used near to the boundaries and at the boundaries are illustrated for the lowest degree case in Fig. 1. Here, the polynomial degree is set to p = 0-and the stencil size inside the domain is n = 2 p + 5 = 5 while the stencil bandwidth can be determined as b = (n − 1)/2. The total number of points of the uniformly structured point-cloud is set to N = 20.

Smooth problem
In the first example, let the distributed line load be chosen as the smooth function f (x) = α(α − 1)x α−2 , where α = 10. Furthermore, the homogeneous BCsũ = 0 m andF = 0 N are given, respectively, at the ends x = 0 and x = L. In this case, the exact solution is a smooth function, i.e. u EX = −x α + α x.
The h-convergence curves obtained for the relative error in 2 -norm of the axial displacement u and force F are plotted against the number of degrees of freedom on log-log scales in Fig. 2. In addition, their rates are also indicated by the slopes. Accordingly, algebraic convergence rates are observed for all refinements. The order of the convergence O(h g ) depends on the polynomial degree p, closely correlated to the theoretical estimate see [4,14,28]. Herein, k denotes the differential order of the considered BVP (k = 2). As expected, higher convergence rates are exhibited for higher-degree polynomial bases. For higher-degrees ( p = 6 and p = 7) however, the stagnation-or, in other words, saturation error-emerges at a certain (but very low) level of the relative error, see [1]. For lower polynomial degrees, the stagnation error is reached at a relatively high error level only with the use of a much larger number of points on the same domain, see [1]. In this example, however, the stagnation error does not yet appear for lower-degree polynomials because the number of points applied is relatively low. Besides, it follows from Eq. (19) that the h-convergence is lost for p = 1, as shown in Fig. 2 as well. The convergence rates are very close to the estimation (19), see also the results for the numerically approximated Laplacian of a smooth function in [14].

Singular problem
In the second example, we consider again a line load f (x) = α(α −1)x α−2 -but now in such a way that the exact solution u EX = −x α + α x is non-smooth. Here, the parameter α controls the smoothness of the solution. If α is small, i.e., α < 1, then the first derivative of the exact solution will become singular at x = 0, see [11]. Thus, we set α to 0.55 in order to study the convergence for a problem with a singular solution. The homogeneous Dirichlet-and Neumann-type BCs,ũ = 0 m andF = 0 N, are specified, respectively, at the boundaries x = 0 and x = L. Again, the h-convergence curves of the relative 2 -error are plotted as the function of the number of degrees of freedom on log-log scales for the displacement u and the axial force F in Fig. 3. The convergence rates are indicated by the corresponding slopes.
The results obtained for p = 1 are not shown since we already observed (in the previous example) that the p-value has to be equal to 2 or higher, as predicted by the theoretical estimate.
In the case of uniform point-cloud refinement, the hconvergence in the displacement error computations is of an algebraic type with the same rate for all polynomial degrees. Concerning the convergence of the axial force, we also observe an algebraic type of convergence with the same rate for different polynomial degrees. The convergence curves are shifted down without changing their slopes as the degree p is increased indicating a lower error constant.
Next, we compare the h-adaptive refinement to the uniform approach. To this end, different polynomial degrees are compared, see Fig. 3. We start the refinement with a uniformly distributed point-cloud with N 0 = 25. It can be observed that the h-adaptive algorithm helps a lot to improve the convergence rates, thereby increasing the accuracy at any considered p-level, not only in the displacement but also in the axial force. A comparison of the numerical and exact solution is given for the 6th iteration step in Fig. 4 for the highest considered polynomial degree p = 5. In this iteration step, the relative 2 -error of the numerical solution computed for the displacement u by the h-adaptive algorithm is already less than 3 %. This error level is much lower than what was achieved by means of the uniform point-distribution.

Shock problem
In the third example, we consider a shock problem with the exact solution which has a high gradient within the domain at x = 1/π , controlled by the parameter α, see [11,32]. Here, α will be set to 120, and f (x) will be defined accordingly. Nonhomogeneous Dirichlet-type BCs are prescribed at both ends of the domain, i.e., at x = 0 and x = L, consistent with the exact solution (20). As in the previous examples, the relative error measured in 2 -norm is plotted in Fig. 5 for both the displacement and the axial force computed with a uniform h-refinement.
In the pre-asymptotic range, it is possible to observe an exponential type of convergence that slows down to an algebraic rate. The convergence rate indicated by the measured slope is almost independent of the chosen polynomial degree.
Henceforward, the adaptive PHS RBF-FD method is studied in 5 refinement steps for p = 2, 3, 4, 5. The results of the adaptive approach are compared to the uniform refinement in Fig. 5. As expected, h-adaptivity significantly improves the convergence rates. The desired accuracy is achieved with much less degrees of freedom for any polynomial degree p, both in the displacement and in the axial force.  To illustrate the numerical approach, the results of the 4th refinement step are once again compared to the uniform approach for the highest investigated polynomial degree p = 5, see Fig. 6. In this adaptive refinement step, the relative 2error of the displacement u and the axial force F is already less than 4 %. Again, the accuracy of the adaptive approach is much higher than that of the uniform refinement.

Two-dimensional boundary value problems
In this section, the convergence behavior of the PHS-based RBF-FD method will be investigated for several representative smooth and non-smooth 2D BVPs of linear elasticity. During the convergence studies, the domain in question is uniformly or quasi-uniformly refined in 6 steps, from 2601 up to 8281 points for all of the benchmark problems. In each step, the degree of the polynomial basis p is kept fixed. In order to observe the influence of the polynomial degree, p will be varied from 2 to 5. Hence, the odd stencil size n has to be adjusted to the actual p-level as n = [( p +1)( p +2)+5] 2 within the domain. The 2D version of the h-adaptive strategy presented in Algorithm 1 will also be tested on the "shock problem" as non-smooth 2D BVP, and its computational performance will be studied.
As an example for the definition of the stencil in the case of uniform point-cloud refinement, a square domain (x, y) ∈ [1, 2] × [1, 2] discretized with a uniformly distributed set of points is considered, see Fig. 7. Different stencils are presented: the symmetric square stencil shape used inside the domain and the half-sided non-symmetric stencil applied next to the boundaries and at the boundaries.
In this example, for the sake of illustration, the stencil size n is set to 5·5 within the domain. The equidistantly structured point set involves N = N I +N B +N E = 15·15+2·15 = 255

Approximation of the two-dimensional differential operators
In the first 2D example, the partial derivatives introduced in Eqs. (2) and (5) are approximated with the RBF-FD method, based on the PHS supplemented with polynomials applying the point-cloud defined on the square domain (x, y) ∈ [1, 2] × [1,2]. Having determined the differential operator matrices D x , D y , and D x x , D yy , as well as D xy , they are applied to the given smooth test function u(x, y) = 1 + sin 4x + cos 3x + sin 2y , Fig. 8 Convergence histories for the relative error in 2 -norm of approximating the first-order differential operators L x and L y acting on the smooth function u Fig. 9 Convergence histories for the relative error in 2 -norm of approximating the second-order differential operators L xx and L yy acting on the smooth function u to compute the derivatives for all points of the structured point-cloud. The 2 -error of the approximation of the firstand second-order partial differential operators, L x , L y , and L x x , L yy , L xy applied to the smooth function (21) is plotted on a double logarithmic scale against the number of degrees of freedom in Figs. 8, 9 and 10. From this, it is evident that the rate of convergence depends on the degree of the polynomial basis p. A higher degree leads to a high algebraic rate of convergence when increasing the number of points. Furthermore, it can be observed that the first-order derivatives can be approximated more accurately as compared to the second-order derivatives.

Smooth problem
Next, the convergence behavior of the PHS RBF-FD method for a unit square plate under linear elastic plane stress conditions is studied. Dirichlet-type BCs are imposed on all of the four boundaries. The body force densities f x , f y and the BCs correspond to the smooth "manufactured" solution: Fig. 10 Convergence histories for the 2 -error-norm of approximating the second-order mixed differential operator L xy acting on the smooth function u u x = − 1 10 cos π (x − 1) sin π (y − 1) , u y = 1 10 sin π (x − 1) 7 sin π (y − 1) 3 . Fig. 11 Convergence curves for the relative error in 2 -norm of the displacements u x and u y -smooth problem Fig. 12 Convergence curves for the relative error in 2 -norm of the stress coordinates σ x , σ y , and τ xy -smooth problem In Figs. 11 and 12 the relative error in 2 -norm of the displacements u x and u y and the stresses σ x , σ y , and τ xy is plotted versus the number of degrees of freedom in double logarithmic style. The stresses are computed by post-processing the first-order partial derivatives of the displacements: where the differential operator matrices D x and D y include the weights for the numerical approximation of the partial derivatives L x and L y at the internal and boundary points, while the column matrices u x and u y consist of the previously-computed displacements u x and u y , respectively. From the computational results, it can be seen that algebraic convergence rates are obtained when increasing the number of points. Again, raising the degree of the additional polynomials improves the convergence rate.
A contour plot of the numerical approximation is compared to the exact solution for the displacements and stresses in Figs. 13 and 14. The numerical solution involves N = 59 · 59 points, a stencil size of n = 17 · 17, and a polynomial degree of p = 2. This comparison demonstrates the good agreement of the numerical approximation with the exact solution.

Shock problem
In our next example, the PHS RBF-FD method will be investigated for the 2D version of the "shock problem" introduced in Sect. 5.1.3. Again, the linear elasticity equation system (2) is numerically solved on a unit square domain, imposing non-homogeneous Dirichlet-type BCs on all parts of the boundary. The body force densities f x , f y and the BCs correspond to the manufactured solution where c = 0.75 m, r 0 = 1 m and α = 60 [11,32].
Since-from an engineering point of view-a sufficiently accurate numerical approximation of the stresses is more of interest than that of the displacements, our numerical investigations are focused on the h-convergence behavior of the stresses. Firstly, in Fig. 15 the relative errors of the uniform hrefinement, measured in 2 -norm, are plotted for the stresses σ x , σ y and τ xy against the number of degrees of freedom in a double logarithmic scale. The relative 2 -error achieved in the 6th refinement step is already below an acceptable level (3 %).
Contour plots of the displacements and stresses are compared to the exact solution in Figs. 16 and 17. The numerical results were obtained with a uniform point-cloud consisting of N = 59 · 59 points, applying a stencil size of n = 17 · 17 and a polynomial degree of p = 2.
Next, the h-adaptive refinement strategy is compared to the uniform approach by extending and applying Algorithm 1 to the numerical solution of 2D BVPs. In the uniform as well as adaptive approach, the refinement procedure is started with a quite coarse resolution, i.e., a uniformly distributed point-cloud containing 7 · 7 = 49 points, keeping the polynomial degree p = 2 and the stencil size n = 17 fixed at each refinement step. However, during the h-adaptive refinement, the stencils become arbitrary-shaped because the knnsearch algorithm is used to find the n-nearest neighbors to the evaluation point of the stencils. Similarly to the 1D BVPs (see Algorithm 1), the h-adaptivity is driven by the following process: The absolute value of the residual is checked for three neighboring points and-when considering the boundaryfor just two neighboring points. A new point is inserted in the center of the neighboring points in question if the refinement criterion (17) along with the definition (18) holds true for one the investigated points.
For quasi-uniform and uniform point-clouds, the convergence behavior of the PHS-RBF methods can be investigated by plotting the relative error curves against either the distance h = 1/ √ N or against N . The convergence rate can be estimated for smooth BVPs by (19). However, for strongly non-uniform point-clouds, the relation between h and N is problem-dependent and the convergence rate can not be determined straightforwardly in each case. Therefore, for the h-adaptive point-cloud refinement strategy and thus also for this 2D comparative study, it is worth measuring the global convergence of relative errors as the function of the so-called effective fill-distance h eff = h /N , where the fill distance     defined as a "mesh size", indicates how well the points in the point-cloud of cover the domain Y ⊆ and represent the radius of the largest empty circle that can be placed among the point locations inside Y [30]. For this reason, the convergence series obtained for the stresses are plotted against the more representative h eff in Fig. 18.
It can also be observed that the 2D version of the hadaptive scheme presented in Algorithm 1 helps a lot to increase both the convergence rates and the accuracy.
To illustrate the h-adaptive refinement approach, the numerical results of the 6th iteration step for the displacements and stresses are plotted on the adaptive point-cloud in Figs. 19 and 20. In this iteration step, the relative 2 -error of the numerical solution computed by the h-adaptive algorithm is already less than 5 %. This error level is lower than what was achieved using the equivalent uniform point-distribution. As a consequence of these results, the extension and the application of this h-adaptive algorithm to the numerical solution of 3D BVPs is the aim of our future research work.

Unstressed circular hole in an infinite plate
Finally, the classical 2D linear elasticity problem of the unstressed circular hole in an infinite plate exposed to unidirectional tension in the x y plane will be studied as a smooth benchmark problem [9,39,40]. Due to symmetry, only one quarter of the plate is investigated. The resulting computational domain is depicted in Fig. 21a.
The exact solution reads for the displacements and σ z = ν σ x + σ y Not only the global but also the local convergence behavior of the RBF-FD method based on the PHS supplemented with a polynomial basis will be tested on this benchmark problem. As a first result, the relative error measured in 2 -norm is plotted for the displacements u x , u y and for the stresses σ x , σ y , τ xy in Figs. 22 and 23. The uniform refinement starts with the point-cloud shown in Fig. 21b. Again, different values for the polynomial degree p are taken into account.
As it is evident from the figures, this setting leads to an algebraic type of convergence that is independent of the degree p.
As a second result, the PHS-based RBF-FD solutions obtained for the stresses σ x , σ y and τ xy at the corners A, B, C, D, and E are depicted together with their exact values are plotted in Figs. 24 and 25. From these figures, quite a fast convergence of pointwise quantities can be observed.
Finally, we present contour plots of the displacements and stresses together with their exact distribution in Figs. 26 and 27. The computations are based on the quasi-uniform point set involving N = 59 · 59 points, applying a stencil size of n = 17 · 17 points together with a polynomial basis of degree p = 2.

Conclusions
A new error-controlled h-adaptive RBF-FD method was developed for the solution of BVPs in linear elasticity. The PHS supplemented with polynomials was chosen as basis.
Firstly, the computational performance of the method was tested through convergence studies of several representative 1D and 2D smooth and non-smooth BVPs. The relative errors were determined in the 2 -norm. It was demonstrated that fast convergence can be obtained by refining the point-cloud, provided that the polynomial degree of the additional polynomials is high enough.
Secondly, the convergence behavior of the newly-developed adaptive point-cloud refinement algorithm was verified in 1D for both the singular and the shock problem, as well as in 2D for the shock problem. The presented h-adaptive scheme leads to a significant improvement both of the convergence rates and the accuracy-not only for the displacement but also for the stress computations. Besides, the advantageous property of the PHS RBF-FD method supplemented with polynomials is that the differentiation weights can be generated in either structured or unstructured (randomized) point-sets.
In the light of the promising numerical results and beneficial properties, the presented h-adaptive point-refinement strategy will be extended to the task of solving 3D BVPs.