Abstract
This paper describes a fully automatic adaptive refinement procedure performed in conjunction with the generalized finite difference method (GFDM) for solving a second-order partial differential equation (PDE) frequently encountered in engineering practice. A quadtree structure is used to organize the clouds of points. The adaptive procedure is based on the existing differences between the gradients of close points. The relative error between two successive adaptation processes is taken as the criterion to stop the algorithm. Using this adaptive procedure, the high gradient regions are easily detected and refined. Some benchmark examples are presented showing the accuracy of the method.
1. Introduction
Considerable research in computational mechanics has been devoted to the development of meshless methods (see Liu [1], Li and Liu [2], and Nguyen et al. [3]). In these methods, the domain of interest is discretized by a scattered set of points. Meshless or meshfree methods can be classified according to the method of solution as weak, local weak, or strong form.
Some meshfree methods have used the weak form. The diffuse element method developed by Nayroles et al. [4] was the first meshless method based on the Galerkin method. Belytschko et al. [5] developed an alternative implementation using moving least squares approximation as defined by Lancaster and Salkauskas [6]. Several other meshless methods such as Partition of Unity Finite Element Method (PUFEM) by Babuška and Melenk [7], h-p cloud by Duarte and Oden [8], and natural element method (NEM) by Sukumar et al. [9] and Cueto et al. [10] have also been published.
Another possibility is to integrate the weak local formulation as in the meshless Petrov-Galerkin method (MLPG), reported by Atluri and Zhu [11, 12].
Another solution is to use the strong form. One of the earliest developments was the SPH method, introduced by Monaghan [13, 14]. Liu et al. [15] proposed a different kind of gridless multiple scale method (RPKM) to improve the accuracy of the SPH method for finite domain problems. Another possibility is to solve the strong form using collocation method as in Oñate et al. [16].
Another important path in the evolution of the strong form has been the development of the generalized finite difference method (GFDM), also called meshless finite difference method. The bases of the GFDM were published by Jensen [17] and Perrone and Kao [18]. Liszka and Orkisz [19] and Orkisz [20] used moving least squares (MLS) in arbitrary irregular grids. Benito et al. [21] developed an implementation of the GFDM showing that the results obtained depend on the number of nodes in the cloud, the relative coordinates of the nodes with respect to the star node, and the weight function employed. Gavete et al. [22] reported improvements of GFDM and compared it with other meshless methods.
Adaptive methods based on the point density increase have been described in Benito et al. [23] and Benito et al. [24, 25] using a posteriori error indicators. A higher order approximation technique based on a modified multipoint approach has also been considered by Jaworska and Orkisz [26].
A totally different adaptive and fully automatic method is proposed in this paper. By combining the GFDM and an adaptive numerical procedure the method can be used to analyze more complicated problems. The method proposed is based on the use of quadtree and on the computations of the gradients, using the differences between the gradients in each quadtree as a simple error indicator. The present approach provides certain restrictions on the location of the nodes due to the use of a quadtree structure. However, the clouds of points generated have an irregular shape. Different partial differential equations of second order are solved in the restrictive case of considering functions with high gradients.
The paper is organized as follows. Section 1 is Introduction. Section 2 describes the GFDM. Section 3 describes the quadtree method used to organize the clouds of points. Section 4 describes the adaptive algorithm and shows the a posteriori error indicator used. Section 5 illustrates the performance of the adaptive method presented for solving second- order different partial differential equations with known analytical solutions. Finally, in Section 6 some conclusions are given.
2. Generalized Finite Difference Method
For any sufficiently differentiable function , in a given domain, the aim is to obtain explicit linear expressions for the approximation of partial derivatives at the points of the domain. First of all, an irregular grid or cloud of points is generated in the domain , where is the boundary of the domain. On defining the central node with a set of nodes surrounding that node, the star then refers to a group of established nodes in relation to a central node. Each node in the domain has an associated star assigned to it. Following [21, 23–25], the explicit finite difference formulae for the first and second spatial derivatives with second-order approximation for the derivatives are obtained: By including the explicit expressions for the values of the partial derivatives , in the equation the star equation is obtained and then with Equation (4) is formed at each point, including the explicit expressions of these derivatives (1) in the Poisson partial differential equation (3). All points in the control scheme are called a star of nodes. The number and the position of nodes in each star are the decisive factors affecting GFD formula approximation. The choice of these supporting nodes is constrained as particular patterns lead to degenerated solutions (Syczewski and Tribillo [27]).
As star selection criterium we follow the denominated cross criterium: the area around the central nodal point, , is divided into four sectors corresponding to quadrants of the Cartesian coordinates system originating at the central node. In each sector two or more nodes are selected, the closest to the origin. If this is not possible, for example, at the boundary, missing nodes can be added to provide the total number of nodes necessary in each star.
Having calculated the values of in the nodes of the domain, we calculate derivatives using formula (1). It is possible to control the precision of GFD solutions by calculating the residual at each point of the interior of the domain using (1) and (2). In order to provide the required and controlled precision of the GFD method, residuals of (2) may be very small and with smoothed distribution over the entire domain. The existence of ill-conditioned stars of nodes depends on the weighting function employed and on the number of nodes by quadrant in each star of nodes, but it can be avoided using the quadtree method shown in the next section.
3. Quadtree Method
A quadtree is a method of organizing data in four quadrants. It is a hierarchical data structure based on the principle of recursive decomposition. A quadtree structure requires each internal node to have exactly four children nodes. When illustrating most quadtree structures, a node that has four child nodes hanging from it is shown, with lines connecting the parent node with its children nodes. The illustration can continue, with four more child nodes hanging from each of the original four child nodes.
The GFDM is a meshfree method; however, it is easy to create an initial simple quadrilateral data structure, as shown in Figure 1(a), as an initial quadtree for the cloud of points. Note that the only restriction used in this paper is that all the quadrilaterals used are convex.
(a)
(b)
An adaptive method using a quadtree can be implemented using an error indicator. Note that the GFDM uses the quadtree structure only to refine the cloud of points; see Figures 1(a) and 1(b) where two quadtree cells are refined. The GFDM with quadtree structure is a fully meshless method. The nodes of each star are selected using the previously defined cross criterion, so there is no relation between the nodes in the quadtree and the nodes in a star of GFDM.
The reason for combining the GFDM and an adaptive numerical procedure is to ameliorate the accuracy of this meshless method for analyzing problems governed by second-order partial differential equations in 2D.
It is interesting to note that not only the number but also the location of the nodes influences the FD operators precision. At the same time it is necessary to take into account that increasing the number of nodes of a cloud of points may decrease the solution error, but it does not change precision of the FD operators.
This organization of a cloud of points using a quadtree structure forces the stars of points used in the GFDM to be in equilibrium; therefore the existence of ill-conditioned stars of nodes is restricted. Also, it all depends on the PDE to be solved.
We can decrease the error by increasing the number of nodes in the cloud. It is well known that the error decreases when the number of nodes is increased in a regular mesh. However, in the case of an irregular mesh, the addition of nodes must be selective, and a solution to this problem is to use the adaptive algorithm described in the following paragraph.
4. Adaptivity Strategy
As is well known, mesh refinement is important in order to capture an accurate solution of a problem. In meshless integral-free methods, a cloud refinement procedure based on gradients can be important and suitable for implementation. Previous gradient based methods have been used by Stein and Rust [28] for finite elements and by Luo and Häussler-Combe [29] for the element-free Galerkin method. The basic idea of the gradient based methods is simple: a larger variation of the gradients needs a richer cloud of nodes.
In this paper a procedure like the one used by Luo and Häussler-Combe [29] has been adapted to the GFDM, but with the following differences:
In each one of the points of the domain, the first derivatives are calculated using (1) and then the gradients at the point are calculated:
In each quadrilateral cell, we calculate the average of the differences of the gradients between the nodes of each quadtree cell (taking into account the distances between the nodes) and we multiply the result by the area of the corresponding cell, to obtain a measure of the gradient intensity in each cell: where
Then we obtain the maximum (max) and minimum (min) values of the cell intensities for the entire domain.
Cells with intensity values greater than (0.75 max intensity + 0.25 min intensity) are selected to be refined, where with and nas = number of adaptation stages. Parameter is related to the number of adaptation stages. Note that when using (9) increases slowly in each refinement step. The quadrilateral cells to be refined are divided into four new cells. A new node is introduced in the midpoint of each nodeless corner of the new cell, respectively (see Figures 1(a) and 1(b), e.g.). The successive clouds of nodes obtained in the refinement process present a smooth nodal density transition.
In the refinement process the algorithm stops when new cells to be refined are not detected. Also the relative error (%) calculated between two successive adaptation processes (10) can be taken as a criterium to stop the algorithm:where is the GFDM solution in node at the adaptive step , is the maximum of the values in the cloud of nodes at the adaptive step , and is the total number of interior nodes of the domain at the adaptive step .
By using this adaptive algorithm in the GFDM, it is possible to decrease the error in the domain by selectively increasing the number of nodes. A fundamental point in using the adaptive procedure is the detection and refinement of the regions with high gradients; to illustrate the procedure some examples have been solved in cases with high gradients in solution.
5. Numerical Results
In this section the previously described h-adaptive method is used to solve various second-order partial differential equations, with the objective of checking the quality of the algorithm. The efficiency of the algorithm has been illustrated analyzing the reduction in the error value of the second-order PDEs with constant or variable coefficients as nodes are added.
The examples have been chosen so that high gradients are present. The potential weighting function used is Parameter has been used in Examples 1–4 and in Example .
As the idea is to analyze the efficiency of the algorithm, it shows that the high gradient regions are well detected and that the adaptive algorithm then adds new nodes selectively, if necessary.
The global exact error is evaluated using the following global error formula: where is the GFDM solution in node at the adaptive step , is the exact value of the solution in the , is the maximum value of the exact values in the cloud of nodes considered, and is the total number of interior nodes of the domain considered.
The following examples show how the h-adaptive method is applied several times, progressively adding new nodes. The nodes added in each step serve as a reference for the following step, in order to obtain more balanced clouds of points at the end of the process. All the examples correspond to different cases with high gradients corresponding to partial differential equations frequently encountered in engineering problems.
5.1. Example 1
For application to solve Poisson equation (13), with Dirichlet boundary conditions being the exact solution (14), can be calculated by substituting exact solution (14) in the partial differential equation (13).
The domain is given in Figure 2 (initial cloud with 50 nodes (32 interior nodes)). After seven refinement stages the final cloud in Figure 3 has 276 nodes (210 interior nodes). The limit used to stop the refinement steps using (10) is 1.5.
In Figure 4 we can see the successive global error calculated using the formula in (12) versus the number of interior nodes.
The error displayed in the logarithmic scale is shown in Figure 5.
5.2. Example 2
The adaptive GFDM is tested in the case of the following partial differential equation: with Dirichlet boundary conditions.
The exact solution is given by The domain is defined by Figure 6 (initial cloud with 25 nodes (9 interior nodes)). Note that in (16) there is a singular point in , not included in the domain. After eight refinement stages the final cloud in Figure 7 has 65 nodes (33 interior nodes). The limit used to stop the refinement steps using (10) has been 0.5.
We can see in Figure 7 how the region with high gradient is detected and refined.
In Figure 8 we can see the successive global errors calculated using formula in (12). The analytical solution corresponds to the case of a function with high gradients in the function and the derivatives.
The error displayed in the logarithmic scale is shown in Figure 9.
5.3. Example 3
The adaptive GFDM is tested with the following partial differential equation with variable coefficients: with Dirichlet boundary conditions.
The exact solution is given by (16) and can be calculated by substituting the exact solution (16) in the partial differential equation (17).
The domain is given in Figure 10 (initial cloud with 43 nodes (17 interior nodes)). Then, after six adaptation stages the refined final cloud with 137 nodes (88 interior nodes) has been obtained; see Figure 11. The limit used to stop the refinement steps using (10) has been 0.75.
In Figure 12 we can see the successive global errors calculated using formula in (12). We can see how the region with high gradient is detected and refined.
The error displayed in the logarithmic scale is shown in Figure 13.
5.4. Example 4
The adaptive GFDM is tested with the following solution of the Poisson equation (18) with Dirichlet boundary conditions.
The exact solution is and can be calculated by substituting the exact solution (18) in the partial differential equation (13).
The domain is given in Figure 14 (initial cloud with 55 nodes (35 interior nodes)). Then, after three adaptation stages, the refined final cloud with 82 nodes (56 interior nodes) is obtained; see Figure 15. The limit used to stop the refinement steps using (10) has been 3.5.
We can see in Figure 15 how the region with high gradients is detected and refined in three adaptation steps.
In Figure 16 we can see the successive global errors calculated using formula in (12). The analytical solution corresponds to the case of a function with high gradients.
The error displayed in the logarithmic scale is shown in Figure 17.
5.5. Example 5
For application to solve a problem taken from the engineering context, a torsion problem defined by the PDE, to calculate the Prandtl stress function in the case of a circular shaft with a circular keyway, see Figure 18.
The exact solution of this problem as given by Saad [30] is
The domain with , is shown in Figure 19 (initial cloud with 86 nodes (60 interior nodes)). After four refinement stages the final cloud of Figure 20 has 1481 nodes (1297 interior nodes). The limit used to stop the refinement steps using (10) is 0.075.
This case is different from the other four above mentioned since the singular strength is lower but it is more widespread in the domain. parameter has been taken as 0.2.
In Figure 21 we can see the successive global error calculated using the formula in (12) versus the number of interior nodes.
The error displayed in the logarithmic scale is shown in Figure 22.
All the examples have been solved using a low initial number of nodes. The h-adaptive method is applied several times progressively adding new nodes, so a smooth transition area is created between zones with different nodal density. As has been shown, the numerical results obtained (see Figures 5, 9, 13, 17, and 22) show a good convergence in a few refinement steps.
6. Conclusions
The use of the GFDM using irregular clouds of points is an interesting way of solving a partial differential equation. This paper shows how a fully automatic adaptive algorithm can be used with the GFDM. The h-adaptive algorithm consists of selectively adding nodes around the areas where the greater errors are located, using a quadtree method to refine the successive clouds of nodes. The algorithm efficiency has been checked analysing the reduction in the solution error value as well as the situation of the added nodes for different partial differential equations, defined in several domains including the cases of a PDE with variable coefficients and a torsion problem.
The computer program developed automatically gives the global estimated error between two successive refinement steps, the number of quadrilaterals with a gradient intensity over the limit, and the number of nodes to be added in each adaptive step. It is then possible to adjust the gradient intensities in the new adaptive steps in order to decrease the error.
Several cases of different second-order partial differential equations, with high gradients in their solutions, have been solved. This shows how the algorithm can detect the high gradient regions and refine them to decrease the error; thus a very good improvement of the accuracy can be obtained by adding a few nodes in each refinement step. All of this can improve the efficiency of the GFDM in future applications.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors acknowledge the support of the Escuela Técnica Superior de Ingenieros Industriales (UNED) of Spain, Project 2014-IFC02.