A Dimension Splitting Generalized Interpolating Element-Free Galerkin Method for the Singularly Perturbed Steady Convection–Diffusion–Reaction Problems

: By introducing the dimension splitting method (DSM) into the generalized element-free Galerkin (GEFG) method, a dimension splitting generalized interpolating element-free Galerkin (DS-GIEFG) method is presented for analyzing the numerical solutions of the singularly perturbed steady convection–diffusion–reaction (CDR) problems. In the DS-GIEFG method, the DSM is used to divide the two-dimensional CDR problem into a series of lower-dimensional problems. The GEFG and the improved interpolated moving least squares (IIMLS) methods are used to obtain the discrete equations on the subdivision plane. Finally, the IIMLS method is applied to assemble the discrete equations of the entire problem. Some examples are solved to verify the effectiveness of the DS-GIEFG method. The numerical results show that the numerical solution converges to the analytical solution with the decrease in node spacing, and the DS-GIEFG method has high computational efﬁciency and accuracy.


Introduction
Since the construction of the approximation function is independent of the mesh, the meshless method, which has developed rapidly in recent years, can completely abandon the mesh reconstruction to ensure high calculation accuracy [1][2][3]. The meshless method has become an essential numerical method in scientific and engineering computation [4][5][6][7].
When studying surface fitting, Lancaster and Salkauskas proposed the moving least squares (MLS) approximation [8], which is based on the traditional least squares (LS) method [9][10][11][12]. The MLS is one of the most important methods for constructing trial functions in meshless methods [13]. Based on the MLS approximation, Belytschko et al. proposed the element-free Galerkin (EFG) method [14], which has a wide range of applications and high calculation accuracy. Since the shape function of the MLS approximation does not satisfy the Kronecker delta property [15], additional numerical techniques are required to impose essential boundary conditions, such as Lagrange multipliers, a penalty method, which may increase the computational burden. For this reason, Lancaster et al. proposed the interpolated moving least-squares (IMLS) method [8] by selecting singular special weight functions. Since the singularity of the weight function is not conducive to numerical calculation, Cheng et al. proposed an improved interpolated moving least squares (IIMLS) method with nonsingular weights for potential problems [16]. The improved meshless method based on the IIMLS has the advantages of nonsingular weight function and direct application of an essential boundary [17,18].

The DS-GIEFG Method for CDR Problems
By using the interpolating shape function of the IIMLS method and coupling the DSM and GEFG method, we will present the DS-GIEFG method for the singularly perturbed steady CDR problems.

The Trial Functions for the DS-GIEFG Method
The trial function of the DS-GIEFG method is presented by using the property of the PU function. From the Refs. [14,16], the following properties are satisfied: where y i is the discrete node, Φ i and b are, respectively, the shape functions and the basis functions used in the IIMLS method [16]. Then, if the basis functions are chosen to be the complete polynomials of order m, the PU function can be established as follows [16] Similar to the GEFG method [22], the trial function of the DS-GIEFG method is defined as where the function τ k (y) is Here, the function θ ij (y) and c ij are, respectively, the nodal enrichment basis function and the corresponding coefficient to be determined of y i , and m i is the number of the local enrichment basis function. The purpose of introducing local enrichment polynomial basis functions is to increase the stability of the meshless method to fluid problems. When m i = 1, the trial function (3) is consistent with that of the EFG method.
According to reference [22,55], we take the local polynomial basis function as the enrichment basis function. In one-dimensional space, θ kj (y) can be taken as the following forms [22,55]: For example, when m i = 2, the trial function of the DS-GIEFG method in onedimensional space is: In Equation (5), the first term represents the effect of the IIMLS method, and the second term is the stabilization parameter of the DS-GIEFG method for the singularly perturbed steady CDR problems.

Discretization Process on the Splitting Plane
The following two-dimensional singular perturbed steady CDR equations are considered: Equation (6) can model the transport of the quantity u = u(x, y). The prescribed function f is the source term describing chemical reactions, and u D is a given functions for the Dirichlet condition. a = (a 1 , a 2 ) T denotes velocity field, 0 <ε 1 denotes the small diffusion coefficient, and c is the reaction coefficient.
Let L denote the number of the splitting plane. According to the idea of the DSM, the problem (6) is divided into a series of one-dimensional problems on the split surface y: where Ω k (k = 1, 2, · · · , L) is a splitting plane such that Define . The weak form of Equation (7) is where , and other symbols are similarly defined, and the inner product is: For Ω (k) , define the numerical solution space as where Φ i (y) is obtained by the IIMLS method on the splitting plane. From Equation (11), the undetermined function u (k) is approximated by where c ij are parameters to be determined. The matrix form of Equation (12) is where c(x) = (c 11 (x), c 12 (x), · · · , c 1m 1 (x), · · · · · · , c n1 (x), c n2 (x), · · · , c nm n (x)) T From Equation (12), the partial derivatives can be obtained as where c ij,x is the partial derivatives concerning x of c (k) ij at x = x k . Then, the discrete Galerkin weak form of Equation (9) is: Substituting Equations (12)- (18) into Equation (19), we can obtain the discrete equations: where c

Assembling the Discrete Equations of the Whole Problem Domain
The discrete equations of all nodes on Ω will be assembled by using the IIMLS method in the x direction. Let where φ i (x) denotes the shape function obtained by the IIMLS method in the x direction with the nodes x k , k = 1, 2, · · · , L, and ∧(x) is an index set of nodes satisfying x in its influence domain. Then, it follows from Equation (25) that Substituting Equations (25)- (27) into Equation (20), it follows On the boundary, the coefficient c ij of enrichment basis function will be supposed to be zero. Then, using the interpolation characteristics of the IIMLS method, and substituting the essential boundary conditions into Equation (28), we can solve the coefficients of all nodes in whole Ω. Additionally, then for ∀(x, y) ∈ Ω, the numerical solution of u is

Numerical Examples
To show the effectiveness of the DS-GIEFG method of this paper, we will present some examples of singularly perturbed steady CDR problems. The examples wer solved on a notebook computer with I5-6200U CPU @2.30GHz, and the SUITESPARSE MATRIX SOFTWARE was also utilized to solve the discrete equations. For comparison, the EFG and VMEFG methods were applied to solve these examples. The weight functions used in this paper were the cubic spline functions, and the influence domains of the nodes in the meshless were rectangular regions with radii of d max × ∆x, where d max is a scalar constant and ∆x denotes the spacing of the nodes. A flow chart for obtaining the numerical solutions of the DS-GIEFG method is listed in Figure 1. The numerical errors are defined as: where u are u h denote the exact and approximated solutions, respectively. Example 1. The first example is a convection-diffusion problem with the following exact solution [56]: The boundary conditions were extracted from the exact solution. When 21 × 21 regular nodes distribution is used, and the parameters are ε = 10 −9 , a 1 = 1 2 , a 2 = −1 3 , c = 0 and d max = 2.2, the numerical solutions solved by the DS-GIEFG method on the lines x = 0.1, 0.2, · · · , 0.9, are listed in Figure 2. We can see that the results of the DS-GIEFG method fit very well with the analytical solutions. For comparison, the EFG and VMEFG methods were also used to solve this example, and the boundary conditions were enforced by the penalty function method with a penalty factor of The influence scalar constants in the EFG and VMEFG methods are Then, under the 21 × 21 regular nodes distribution, the absolute errors on the whole field are shown in Figure 3. Clearly, the solution of the DS-GIEFG method has better calculation accuracy than those of the other two methods. For comparison, the EFG and VMEFG methods were also used to solve this example, and the boundary conditions were enforced by the penalty function method with a penalty factor of α = 10 5 . The influence scalar constants in the EFG and VMEFG methods are d max = 2.8. Then, under the 21 × 21 regular nodes distribution, the absolute errors on the whole field are shown in Figure 3. Clearly, the solution of the DS-GIEFG method has better calculation accuracy than those of the other two methods.
To show the convergence on node spacing, when 11 × 11, 21 × 21, · · · , 51 × 51 regular node distributions were used; the relative errors of the DS-GIEFG, VMEFG and EFG methods with respect to node spacing are shown in Figure 4, and the corresponding CPU times are listed in Figure 5. Figure 4 shows that the proposed method has higher calculation accuracy and convergence order, and the numerical solution will gradually converge to the analytical solution with the reduction in node spacing. It can also be seen from Figure 5 that the proposed method has less computation time and higher computational efficiency than the EFG and VMEFG methods. For comparison, the EFG and VMEFG methods were also used to solve this example, and the boundary conditions were enforced by the penalty function method with a penalty factor of 5 10 α = . The influence scalar constants in the EFG and VMEFG methods are Then, under the 21 × 21 regular nodes distribution, the absolute errors on the whole field are shown in Figure 3. Clearly, the solution of the DS-GIEFG method has better calculation accuracy than those of the other two methods.  To show the convergence on node spacing, when 11 × 11, 21 × 21, ,51 × 51 regular node distributions were used; the relative errors of the DS-GIEFG, VMEFG and EFG methods with respect to node spacing are shown in Figure 4, and the corresponding CPU times are listed in Figure 5. Figure 4 shows that the proposed method has higher calculation accuracy and convergence order, and the numerical solution will gradually converge to the analytical solution with the reduction in node spacing. It can also be seen from Figure  5 that the proposed method has less computation time and higher computational efficiency than the EFG and VMEFG methods.  thematics 2021, 9, 2524 9 (c) EFG To show the convergence on node spacing, when 11 × 11, 21 × 21, ,51 × 51 reg node distributions were used; the relative errors of the DS-GIEFG, VMEFG and EFG m ods with respect to node spacing are shown in Figure 4, and the corresponding CPU ti are listed in Figure 5. Figure 4 shows that the proposed method has higher calcula accuracy and convergence order, and the numerical solution will gradually converg the analytical solution with the reduction in node spacing. It can also be seen from Fig  5 that the proposed method has less computation time and higher computational ciency than the EFG and VMEFG methods.     Figure 7. Additionally, the absolute er of the DS-GIEFG and VMEFG methods on the whole field are shown in Figure 8. It can seen from Figure 6a that the solution of the EFG method is unstable, and there is a n physical oscillation for the diffusion convection problem with a very small diffusion c ficient. Figure 6 shows that both VMEFG and DS-GIEFG methods can obtain stable merical solutions for small diffusion coefficients. Figure 8 also shows that the DS-GIE method can have higher calculation accuracy than the VMEFG method.

Example 2.
This example [57] is defined in the field of Ω = [0, 1] × [0, 1] with the exact solution of: The essential boundary was used and determined by the exact solution.
The parameters were chosen to be ε = 10 −8 , a 1 = a 2 = 1, c = 0 and d max = 2.2. When using 21 × 21 node distributions, the 3D surfaces of the solution obtained by the EFG, VMEFG and DS-GIEFG methods are, respectively, shown in Figure 6a-c, and the corresponding exact solution is plotted in Figure 6d. The penalty factor and scalar constant of influence domain in the EFG and VMEFG methods are α = 10 5 and d max = 2.6, respectively. The corresponding snapshots of the solutions of the DS-GIEFG method on the fields of lines x = 0.1, 0.2, · · · , 0.9 are listed in Figure 7. Additionally, the absolute errors of the DS-GIEFG and VMEFG methods on the whole field are shown in Figure 8. It can be seen from Figure 6a that the solution of the EFG method is unstable, and there is a non-physical oscillation for the diffusion convection problem with a very small diffusion coefficient. Figure 6 shows that both VMEFG and DS-GIEFG methods can obtain stable numerical solutions for small diffusion coefficients. Figure 8 also shows that the DS-GIEFG method can have higher calculation accuracy than the VMEFG method.

Example 2.
This example [57] is defined in the field of The essential boundary was used and determined by the exact solution.
The parameters were chosen to be  Figure 7. Additionally, the absolute errors of the DS-GIEFG and VMEFG methods on the whole field are shown in Figure 8. It can be seen from Figure 6a that the solution of the EFG method is unstable, and there is a nonphysical oscillation for the diffusion convection problem with a very small diffusion coefficient. Figure 6 shows that both VMEFG and DS-GIEFG methods can obtain stable numerical solutions for small diffusion coefficients. Figure 8 also shows that the DS-GIEFG method can have higher calculation accuracy than the VMEFG method.  When the 17 × 17, 33 × 33, 65 × 65, 129 × 129 node distributions are used, the relative errors of the EFG, VMEFG and DS-GIEFG methods in the lg-lg scale are presented in Figure 9. The results of Ref. [57] are also listed in this figure. Additionally, the corresponding CPU time of these methods is plotted in Figure 10. For this example, there is no doubt that the DS-GIEFG method in this paper has higher calculation accuracy and efficiency.  When the 17 × 17, 33 × 33, 65 × 65, 129 × 129 node distributions are used, the relative errors of the EFG, VMEFG and DS-GIEFG methods in the lg-lg scale are presented in Figure 9. The results of Ref. [57] are also listed in this figure. Additionally, the corresponding CPU time of these methods is plotted in Figure 10. For this example, there is no doubt that the DS-GIEFG method in this paper has higher calculation accuracy and efficiency.  When the 17 × 17, 33 × 33, 65 × 65, 129 × 129 node distributions are used, the relative errors of the EFG, VMEFG and DS-GIEFG methods in the lg-lg scale are presented in Figure 9. The results of Ref. [57] are also listed in this figure. Additionally, the corresponding CPU time of these methods is plotted in Figure 10. For this example, there is no doubt that the DS-GIEFG method in this paper has higher calculation accuracy and efficiency. When the 17 × 17, 33 × 33, 65 × 65, 129 × 129 node distributions are used, the relative errors of the EFG, VMEFG and DS-GIEFG methods in the lg-lg scale are presented in Figure 9. The results of Ref. [57] are also listed in this figure. Additionally, the corresponding CPU time of these methods is plotted in Figure 10. For this example, there is no doubt that the DS-GIEFG method in this paper has higher calculation accuracy and efficiency. Ref. [57] DS-GIEFG EFG VMEFG Figure 9. The variation of relative errors of different methods with respect to the node spacin lg-lg scale).   Ref. [57] DS-GIEFG EFG VMEFG Figure 9. The variation of relative errors of different methods with respect to the node spacing lg-lg scale).

Example 3.
This example is chosen to be an interior layer problem with the exact solution [58]: Let ε = 10 −9 , a 1 = 0, a 2 = 1, c = 1, a = b = 0.05. The other parameters are exactly the same as in Example 2. When using 51 × 51 nodes, the 3D surface of the exact and numerical solutions of the DS-GIEFG method are plotted in Figure 11. Additionally, the snapshots of solutions on the lines y = 0.1, 0.2, · · · , 0.5 are shown in Figure 12. These two figures show that the solutions of the DS-GIEFG method are very consistent with those of the exact solutions.
When the 9 × 9, 17 × 17, 33 × 33, 65 × 65, 129 × 129 nodes distributions are used, the errors u − u h of the Ref. [58], EFG, VMEFG and DS-GIEFG methods are listed in Table 1. Additionally, the corresponding variation of the errors with respect to the CPU time is shown in Figure 13. These results once again show that the method proposed in this paper has high computational accuracy and computational efficiency.  Table 1. Additionally, the corresponding variation of the errors with respect to the CPU time is shown in Figure 13. These results once again show that the method proposed in this paper has high computational accuracy and computational efficiency.  When the 9×9, 17×17, 33×33, 65×65, 129×129 nodes distributions are used, the er h || || u u − of the Ref. [58], EFG, VMEFG and DS-GIEFG methods are listed in Tab Additionally, the corresponding variation of the errors with respect to the CPU tim shown in Figure 13. These results once again show that the method proposed in this pa has high computational accuracy and computational efficiency.

Conclusions
To improve the computational efficiency and overcome the non-physical oscilla of the EFG method in solving singularly perturbed fluid problems, a DS-GIEFG met for singularly perturbed steady CDR problems is proposed in this paper by construc the trial functions based on the IIMLS method and coupling the DSM and GEFG meth The stabilization parameter of the DS-GIEFG method is only related to the nodes and

Conclusions
To improve the computational efficiency and overcome the non-physical oscillation of the EFG method in solving singularly perturbed fluid problems, a DS-GIEFG method for singularly perturbed steady CDR problems is proposed in this paper by constructing the trial functions based on the IIMLS method and coupling the DSM and GEFG method. The stabilization parameter of the DS-GIEFG method is only related to the nodes and not to the mesh. Since the two-dimensional problem is divided into a series of one-dimensional problems, the DS-GIEFG method has high efficiency. Three numerical examples are solved by the DS-IEFG method of this paper. The numerical solutions are compared with those of the EFG and VMEFG methods. The numerical results show that the solutions of the DS-GIEFG method converge to the analytical solution with the increase in the number of nodes, and the DS-GIEFG method has higher computational efficiency and accuracy than the EFG and VMEFG methods.
Author Contributions: Conceptualization, F.S. and J.W.; methodology, F.S. and J.W.; software, J.W. and X.K.; writing-original draft preparation, J.W. and F.S.; writing-review and editing, F.S. and R.C. All authors have read and agreed to the published version of the manuscript.