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

Abstract

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 efficiency and accuracy.

1. 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].

Many meshless methods have been presented [19], such as the smooth particle hydrodynamics method, element-free Galerkin (EFG) method, reproducing kernel particle method, hp-cloud method, and meshless local Petrov–Galerkin (MLPG) method. Peridynamics (PD) is also an important meshless method for solving solid and fluid mechanics problems [20,21]. When the EFG method is used to solve some fluid problems, the numerical solution is prone to non-physical oscillation [22,23,24]. By coupling the variational multiscale method and EFG method, Ouyang et al. established the variational multiscale EFG method (VMEFG) [25]. The VMEFG method has been successfully applied in many physical problems [26,27,28]. However, the disadvantage of VMEFG is that the bubble function in the stabilization factor must be based on the mesh and is difficult to construct in high-dimensional cases. Li et al. proposed the general element-free Garlerkin (GEFG) method [22] for Stokes problems. The GEFG method overcomes the non-physical oscillation of the solution caused by the EFG method, and also overcomes the deficiency that the stabilization factor of the VMEFG method must be related to the mesh.

Due to the complexity of the shape function, the traditional EFG method requires a large amount of calculation, especially for three-dimensional problems. Then, Cheng et al. proposed the dimension-splitting element-free Galerkin method by coupling the dimension-splitting method (DSM) and EFG method [29,30,31]. Since the DSM divides high-dimensional problems into a series of low-dimensional problems and then iteratively solves them, it can effectively improve the computational efficiency of numerical methods [32,33]. Therefore, the dimension-splitting element-free Galerkin method can significantly improve the computational efficiency of the EFG method [34].

The convection–diffusion–reaction (CDR) equation has a wide range of applications in the fields of physics and chemistry [35,36,37]. The steady CDR equation can model the transport of the quantity and represent the balance of three processes: transport by convection, diffusion and reactions [37]. As it is difficult to obtain the analytical solutions of the CDR equation in complex situations, many numerical methods have been proposed to obtain the numerical solutions [38,39,40]. A finite difference method [36] and two spurious oscillations at layers diminishing (SOLD) methods [37] have been presented for solving steady CDR problems. For the singularly perturbed problems, the numerical solutions are prone to instability issues [41,42]. Then, when the diffusion coefficient is very small in CDR problems, it is usually necessary to add additional stability measures to ensure that the numerical method obtains a stable solution [43]. Wu et al. presented an adaptive upwind finite difference method for the CDR problem with two small parameters [44]. The finite difference method [45], a cubic B-spline-based semi-analytical algorithm [46], and weak Galerkin method [47] are proposed for studying the numerical solutions of the CDR problems.

Due to the advantages of the meshless method, some meshless methods have been developed for the CDR problems [48,49]. Zhang and Xiang discussed the numerical solutions of the CDR equation with small diffusion by the VMEFG method [50,51]. Li et al. also applied the local Petrov Galerkin method to the multi-dimensional CDR equations based on the radial basis function [52]. Additionally, the local knot method [53] and the meshless boundary collocation method [54] have been presented to solve the numerical solution of the CDR equation.

In this paper, by constructing the trial functions from the IIMLS method and coupling the DSM and GEFG method, a dimension splitting general interpolating element-free Galerkin (DS-GIEFG) method is presented for the singularly perturbed steady CDR problems. In the DS-GIEFG method, DSM divides the two-dimensional problem into a series of one-dimensional problems. The GEFG method is used to establish the discrete equations on the low-dimensional problems. Finally, the IIMLS method is used to assemble the discrete equations on the entire problem domain. The DS-GIEFG method will have high computational efficiency and solution stability for the singularly perturbed steady CDR problem.

2. 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.

2.1. 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:

$${\displaystyle \sum _{i=1}^n{\Phi }_i(y)b(y_i)}=b(y)$$

(1)

where $y_i$ is the discrete node, ${\Phi }_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]

$${\displaystyle \sum _{i=1}^n{\Phi }_i}=1$$

(2)

Similar to the GEFG method [22], the trial function of the DS-GIEFG method is defined as

$$u^h(y)={\displaystyle \sum _{i=1}^n{\Phi }_i(y){\tau }_i(y)}$$

(3)

where the function ${\tau }_k(y)$ is

$${\tau }_i(y)={\displaystyle \sum _{j=1}^{m_i}{\theta }_{ij}(y)c_{ij}}$$

(4)

Here, the function ${\theta }_{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, ${\theta }_{kj}(y)$ can be taken as the following forms [22,55]:

(a)

If $m_i=1$, then ${\theta }_{i1}(y)=1$.

(b)

if $m_i=2$, then $\left\{{\theta }_{i1},{\theta }_{i2}\right\}=\left\{1,{(y-y_i)}^2\right\}$.

(c)

if $m_i=3$, then $\left\{{\theta }_{i1},{\theta }_{i2},{\theta }_{i3}\right\}=\left\{1,{(y-y_i)}^2,{(y-y_i)}^3\right\}.$

For example, when $m_i$ = 2, the trial function of the DS-GIEFG method in one-dimensional space is:

$$u^h={\displaystyle \sum _{i=1}^n{\Phi }_i(y)\cdot c_{i1}}+{\displaystyle \sum _{i=1}^n{\Phi }_i(y)\cdot {(y-y_i)}^2{c}_{i2}}$$

(5)

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.

2.2. Discretization Process on the Splitting Plane

The following two-dimensional singular perturbed steady CDR equations are considered:

$$\{\begin{array}{ll}a\cdot \nabla u-\epsilon \mathsf{\Delta }u+cu=f,\hfill & \mathrm{in}\Omega \subset R^2,\hfill \\ u=u_D,\hfill & \mathrm{on}\mathrm{boundary}\Gamma .\hfill \end{array}$$

(6)

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)}^{\mathrm{T}}$ denotes velocity field, $0<\epsilon \ll 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$:

$$\{\begin{array}{ll}{\left(a\cdot \nabla u-\epsilon \mathsf{\Delta }u+cu\right)|}_{x=x_k}={f|}_{x=x_k},\hfill & \mathrm{in}{\Omega }_k=\Omega \cap \left\{x=x_k\right\},\hfill \\ u=u_D,\hfill & \mathrm{on}\mathrm{boundary}\Gamma ,\hfill \end{array}$$

(7)

where ${\Omega }_k$ ($k=1,2,\cdots ,L$) is a splitting plane such that

$$\Omega =\underset{k=1}{\overset{L}{\cup }}{\Omega }_k\times (x_{k-1},x_k]$$

(8)

Define $V^{(k)}\triangleq H^1({\Omega }_k)\cap C^0({\Omega }_k)$. The weak form of Equation (7) is

$$\left(a\cdot \nabla u^{(k)},v\right)+\epsilon \left(v_y,u_y^{(k)}\right)-\epsilon \left(v,u_{xx}^{(k)}\right)+c\left(v,u^{(k)}\right)=(v,f^{(k)}),\forall v\in V^{(k)}$$

(9)

where $u^{(k)}=u{|}_{x=x_k}$,$u_y^{(k)}=\frac{\partial u}{\partial y}{|}_{x=x_k}$, $u_{xy}^{(k)}=\frac{{\partial }^{2}u}{\partial x\partial y}{|}_{x=x_k}$, and other symbols are similarly defined, and the inner product is:

$$(v,u)\triangleq {\displaystyle {\int }_{{\Omega }_k}uvdy}$$

(10)

For ${\Omega }^{(k)}$, define the numerical solution space as

$$V^h=\mathrm{span}\{{\Phi }_i(y){\theta }_{ij}(y),1\le i\le n,1\le j\le m_i,y\in {\Omega }^{(k)}\}$$

(11)

where ${\Phi }_i(y)$ is obtained by the IIMLS method on the splitting plane.

From Equation (11), the undetermined function $u^{(k)}$ is approximated by

$$u^{(k)}\approx {\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{m_i}{\Phi }_i(y){\theta }_{ij}(y)c_{ij}^{(k)}}}$$

(12)

where $c_{ij}^{(k)}$ are parameters to be determined.

The matrix form of Equation (12) is

$$u^{(k)}={\varpi }^{(k)}{c}^{(k)}(x)$$

(13)

where

$${\varpi }^{(k)}={({\Phi }_1(y){\theta }_1^{(k)}(y),{\Phi }_2(y){\theta }_2^{(k)}(y),\cdots ,{\Phi }_n(y){\theta }_n^{(k)}(y))}^{\mathrm{T}}$$

(14)

$${\theta }_i^{(k)}=({\theta }_{i1}(y),{\theta }_{i2}(y),\cdots ,{\theta }_{i{m}_i}(y))$$

(15)

$$c(x)={(c_{11}(x),c_{12}(x),\cdots ,c_{1m_1}(x),\cdots \cdots ,c_{n1}(x),c_{n2}(x),\cdots ,c_{n{m}_n}(x))}^{\mathrm{T}}$$

(16)

From Equation (12), the partial derivatives can be obtained as

$$u_y^{(k)}\approx {\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{m_i}\frac{\partial }{\partial y}[{\Phi }_i(y){\theta }_{ij}(y)]c_{ij}^{(k)}}}$$

(17)

$$u_x^{(k)}\approx {\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{m_i}{\Phi }_i(y){\theta }_{ij}(y)c_{ij,x}^{(k)}}}$$

(18)

where $c_{ij,x}^{(k)}$ is the partial derivatives concerning $x$ of $c_{ij}^{(k)}$ at $x=x_k$.

Then, the discrete Galerkin weak form of Equation (9) is: for $\forall {\overline{u}}^{(k)}\in V^h,$ find ${\widehat{u}}^{(k)}\in V^h$ such that

$$\left({\overline{u}}^{(k)},a\cdot \nabla {\widehat{u}}^{(k)}\right)+\epsilon \left({\overline{u}}_y^{(k)},{\widehat{u}}_y^{(k)}\right)-\epsilon \left({\overline{u}}^{(k)},{\widehat{u}}_{xx}^{(k)}\right)+c\left({\overline{u}}^{(k)},{\widehat{u}}^{(k)}\right)=({\overline{u}}^{(k)},f^{(k)})$$

(19)

Substituting Equations (12)–(18) into Equation (19), we can obtain the discrete equations:

$$-\epsilon M_0^{(k)}{c}_{xx}^{(k)}+a_1{M}_0^{(k)}{c}_x^{(k)}+\left[\epsilon M_2^{(k)}+c{M}_0^{(k)}+a_2{M}_1^{(k)}\right]c^{(k)}=F^{(k)}$$

(20)

where $c_{xx}^{(k)}=\frac{{\partial }^{2}c(x)}{\partial x^2}{|}_{x=x_k}$, $c_x^{(k)}=\frac{\partial c(x)}{\partial x}{|}_{x=x_k}$,$c^{(k)}=c(x_k)$, and

$$M_2^{(k)}={\displaystyle {\int }_{{\Omega }^{(k)}}\frac{\partial }{\partial y}{\varpi }^{(k)}\cdot {\left(\frac{\partial }{\partial y}{\varpi }^{(k)}\right)}^{\mathrm{T}}dy}$$

(21)

$$M_1^{(k)}={\displaystyle {\int }_{{\Gamma }^{(k)}}{\varpi }^{(k)}\cdot {\left(\frac{\partial }{\partial y}{\varpi }^{(k)}\right)}^{\mathrm{T}}dy}$$

(22)

$$M_0^{(k)}={\displaystyle {\int }_{{\Omega }^{(k)}}{\varpi }^{(k)}\cdot {\left({\varpi }^{(k)}\right)}^{\mathrm{T}}dy}$$

(23)

$$F^{(k)}={\displaystyle {\int }_{{\Omega }^{(k)}}{\varpi }^{(k)}\cdot f^{(k)}dy}$$

(24)

2.3. Assembling the Discrete Equations of the Whole Problem Domain

The discrete equations of all nodes on $\Omega$ will be assembled by using the IIMLS method in the $x$ direction. Let

$$c(x)={\displaystyle \sum _{s\in \wedge (x)}{\varphi }_s(x)c(x_s)}={\displaystyle \sum _{s\in \wedge (x)}{\varphi }_s(x)c^{(s)}}$$

(25)

where ${\varphi }_i(x)$ denotes the shape function obtained by the IIMLS method in the $x$ direction with the nodes $x_k$, $k=1,2,\cdots ,L$, and $\wedge (x)$ is an index set of nodes satisfying $x$ in its influence domain.

Then, it follows from Equation (25) that

$$c_x^{(k)}={\displaystyle \sum _{s\in \wedge (x_k)}\frac{\partial }{\partial x}{\varphi }_s(x_k)c^{(s)}}$$

(26)

$$c_{xx}^{(k)}={\displaystyle \sum _{s\in \wedge (x_k)}\frac{{\partial }^2}{\partial x^2}{\varphi }_s(x_k)c^{(s)}}$$

(27)

Substituting Equations (25)–(27) into Equation (20), it follows

$$\{\begin{array}{l}\left[\epsilon M_2^{(k)}+c{M}_0^{(k)}+a_2{M}_1^{(k)}\right]c^{(k)}+{\displaystyle \sum _{s\in \wedge (x_k)}\left(a_1\frac{\partial }{\partial x}{\varphi }_s(x_k)M_0^{(k)}-\epsilon \frac{{\partial }^2}{\partial x^2}{\varphi }_s(x_k)M_0^{(k)}\right)}{c}^{(s)}=F^{(k)}\hfill \\ k=1,2,\cdots ,L.\hfill \end{array}$$

(28)

On the boundary, the coefficient $c_{ij}^{(k)}$ 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 $\Omega$. Additionally, then for $\forall (x,y)\in \Omega$, the numerical solution of $u$ is

$$u^h(x,y)={\displaystyle \sum _{s\in \wedge (x)}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{m_i}{\varphi }_s(x){\Phi }_i(y){\theta }_{ij}(y)c_{ij}^{(s)}}}}$$

(29)

3. 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 }\times \mathsf{\Delta }x$, where $d_{\max }$ is a scalar constant and $\mathsf{\Delta }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:

$$\left|\right|\mathrm{RE}\left|\right|=\frac{\left|\right|u-u^{\mathrm{h}}|{|}_{L^2(\Omega )}}{\left|\right|u|{|}_{L^2(\Omega )}}$$

(30)

$$\left|\right|e\left|\right|=\left|\right|u-u^{\mathrm{h}}|{|}_{L^2(\Omega )}$$

(31)

where $u$ are $u^{\mathrm{h}}$ denote the exact and approximated solutions, respectively.

Example 1.

The first example is a convection–diffusion problem with the following exact solution [56]:

$$u(x,y)=\frac{\mathrm{65,536}}{729}{x}^3(1-x)y^3(1-y), (x,y)\in \Omega =[0,1]\times [0,1]$$

(32)

The boundary conditions were extracted from the exact solution.

When 21 × 21 regular nodes distribution is used, and the parameters are $\epsilon ={10}^{-9}$, $a_1=\frac{1}{2},a_2=\frac{-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, $\cdots$, 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 $\alpha ={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, $\cdots$, 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.

Example 2.

This example [57] is defined in the field of $\Omega =[0,1]\times [0,1]$ with the exact solution of:

$$u(x,y)=100x^{2}y{(1-x)}^2(1-y)(1-2y)$$

(33)

The essential boundary was used and determined by the exact solution.

The parameters were chosen to be $\epsilon ={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 $\alpha ={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, $\cdots$, 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.

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.

Example 3.

This example is chosen to be an interior layer problem with the exact solution [58]:

$$u(x,y)=\frac{1}{2}x(1-x)y(1-y)(1-\tan \frac{y-a}{b}),(x,y)\in \Omega =[0,1]\times [0,1]$$

(34)

Let $\epsilon ={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,\cdots ,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 $\left|\right|u-u^{\mathrm{h}}\left|\right|$ of the Ref. [58], EFG, VMEFG and DS-GIEFG methods are listed in

Table 1. The errors $\left|\right|u-u^{\mathrm{h}}\left|\right|$ of the Ref. [58], EFG, VMEFG and DS-GIEFG methods for different nodes spacing.

Nodes	DS-GIEFG	VMEFG	EFG	Ref. [58]
9 × 9	7.90 × 10−4	1.10 × 10−1	4.02 × 10−1	3.80 × 10−3
17 × 17	1.81 × 10−4	1.36 × 10−2	1.14 × 10−1	1.31 × 10−3
33 × 33	1.37 × 10−5	4.74 × 10−4	5.29 × 10−3	4.37 × 10−4
65 × 65	9.46 × 10−8	6.91 × 10−5	3.89 × 10−4	1.16 × 10−4
129 × 129	1.29 × 10−8	1.52 × 10−5	7.75 × 10−5	2.96 × 10−5

. 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.

4. 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.