An Adaptive Domain Partitioning Technique for Meshfree-Type Methods

Abstract

An overlapping domain partitioning based on adapting nodes is presented for the meshless-type methods. The decomposition of the domain is carried out based on the distribution of the nodes produced rather than the geometry of the problem. A set of adaptive nodes is first generated using the dimension reduction and equidistributing along the coordinate directions with respect to arc-length monitor. The domain is then partitioned in such a way that the same number of nodes are allocated to the subdomains. A radial basis function collocation method is applied to each subdomain followed by assembling the global solution from the subproblem's solutions. A generalized thin plate spline with sufficient smoothness is used as a basis function in the collocation method. Some numerical results will be presented to show the performance of the proposed method.

1. Introduction

During the last two decades, numerical solutions of partial differential equations (PDEs) based on radial basis functions (RBFs) have received a considerable attention. This method was first introduced by Kansa [1] as an unsymmetric collocation method followed by several approaches such as the method of fundamental solutions [2] and the Galerkin’s method [3]. These methods are known as mesh-free methods due to their advantage of avoiding mesh construction. The meshless methods are often divided into two major categories: the boundary-type method (see e.g., [2, 4]) and the domain-type method (see e.g., [5, 6]). The current paper aims to improve the Kansa’s method which belongs to the latter category, though the idea can be applied to the former category as well.

The RBFs are a set of powerful basis functions for multivariate approximation with many advantages including high accuracy, ease of implementation, and low memory storage. However, they suffer from ill-conditioning, specially, when the size of the interpolation matrix increases. Moreover, the method involves dense matrices, which reduces the computational efficiency of the underlying method. To avoid the above difficulties, some techniques such as domain decomposition methods (DDM) [5, 7], compactly supported RBFs (CS-RBFs) [8], and preconditioning [9] have been recommended. In addition, a mesh-free numerical procedure based on the local radial basis function collocation method with an adaptive nodes technique was presented in [10].

To improve the conditioning of the interpolation matrix, we suggest an overlapping DDM based on an adaptive nodes strategy rather than only the geometry of the problem. More precisely, a set of adaptive nodes is first produced in the domain and then the domain is partitioned into some overlapping subregions such that the same number of nodes are located in the subdomains. The classical DDM Schwarz method is applied to the partitioned problem, that is, applying the collocation meshless method to each subproblem and obtaining the local solutions followed by updating the interface boundary solutions. This process is performed iteratively until a desired accuracy for the global solution is achieved. Using this method, we take advantages of both DDM and adaptive nodes method. On one hand, the method takes benefit of the adaptive mesh method which aims to use a minimum number of nodes. On the other hand, it enjoys the advantages of the DDM, which reduce the size of the matrices involved and, consequently, improve the conditioning of the problem.

As noted, an adaptive nodes technique is required for the proposed method. It is well known that the main idea in adapting meshes is to use a minimum number of nodes while keeping a desired accuracy. This is achieved by allocating more mesh points to the areas where the solution has rapid variation. There have been many adaptive mesh strategies in the literature [11, 12]. We employ an equidistribution strategy based on the dimension reduction which permits mesh points generation in the higher dimensions [13]. In particular, we use a method presented in [14] for rectangular domains based on dimension reduction and equidistributing along the grid lines in the coordinate directions. In the dimension reduction method, the equidistribution process is reduced to a 1D case. Also a generalization of the above-mentioned method to the case of three dimensions has been presented in [15]. In addition, an extension of the 2D adaptive mesh to the case of irregular domains has been presented in [16].

As is well known, in the equidistribution strategy, the mesh distribution is carried out in such a way that some measure of error, called a monitor function, is equalized over each subinterval. In the current work, an adaptive nodes technique is used to find a domain partitioning in such a way that the solution error is roughly equalized over each local problem. This paper is organized as follows. In Section 2, an adaptive mesh technique used in this study is introduced. The overlapping Schwarz method will be briefly described in Section 3. The new adaptive partitioning method is presented in Section 4. In Section 5, the collocation meshless method is reviewed. Some numerical results will be presented in Section 6.

2. Adaptive Nodes in a Rectangle

For a 2D domain Ω, a partition {Ω_𝑖}^𝑛𝑖=1 is called equidistributing with respect to a nonnegative piecewise continuous function 𝑀(𝑥,𝑦) (see [17]) if_Ω𝑖𝑀(𝑥,𝑦)𝑑𝑥𝑑𝑦=constant,(2.1) where the monitor function 𝑀 depends on the solution of the underlying PDE and its derivatives. To construct suitable monitor functions for different applications is a line of research [18]; however, the current work focuses on the mesh generation strategy using a well-known monitor, arc-length, which is given by (see [19])𝑀=1+𝑢^2𝑥.(2.2)

We now briefly review a 2D adaptive mesh method presented in [14] based on dimension reduction and equi-distributing along the coordinate axes. The method starts with a uniform mesh in a rectangle in the form{(𝑥,𝑦)∣𝑎≤𝑥≤𝑏,𝑐≤𝑦≤𝑑}.(2.3) The equidistribution process is performed in three stages. In the first stage, equidistribution is performed in the direction of the 𝑥-axis (Figure 3(a)). In the second stage, the equidistribution is performed in the vertical direction along the grid lines produced in the first stage. Since the grid lines are curved, the distribution is performed along the arc rather than the vertical coordinate direction (Figure 3(b)). A similar procedure is performed in the third stage along the horizontal grid curves, using the monitor in the 𝑥 coordinate direction again (Figure 1(c)).

(a)

(b)

(c)

Figure 1 

The three stages of the adaptive mesh method are shown in Figures (a), (b), and (c), respectively, and in each direction two grid curves are displayed with the quadrilateral formed.

The mesh resulting from the above procedure forms quadrilaterals whose sides equi-distribute the grid lines in the two coordinate directions (see Figure 1(c)). The produced mesh may not be precisely equi-distributing in 2D; however, the mesh points are suitable for any meshless method in which the mesh points, rather than mesh lines, are important. It should be noted that the connectivity of the mesh, which is used in the mesh generation process, is not utilized in our proposed adaptive ddm, which is applied to a meshless method.

This mesh generation technique has two important properties, which will be used in the current work.(1)It keeps the number of initial uniform mesh points by which the adapting process is performed. More precisely, if 𝐴=𝑥_𝑖,𝑦_𝑖∶𝑖=1,…,𝑁(2.4) is the initial set of equally spaced nodes and 𝐵=𝑥^𝑖,𝑦^𝑖∶𝑖=1,…,𝑁^(2.5) is the set of points produced by the adapting strategy, then 𝑁=𝑁^ and there will be a one-to-one correspondence between the members of 𝐴 and 𝐵. In fact, the adaptive method moves a point (𝑥_𝑖,𝑦_𝑖) to a new position denoted by (𝑥^𝑖,𝑦^𝑖).(2)It keeps the order of points, that is, for every two points (𝑥_𝑖,𝑦_𝑖) and (𝑥_𝑗,𝑦_𝑗), if 𝑥_𝑖<𝑥_𝑗, then 𝑥^𝑖<𝑥^𝑗, and if 𝑦_𝑖<𝑦_𝑗, then 𝑦^𝑖<𝑦^𝑗. The 2D adaptive method inherits these properties from the equi-distributing algorithm in the case of 1D presented in [17] taking into account that the dimension reduction is used to reduce the problem to the 1D case.

This technique will be used in Section 4 to combine the DDM and the adaptive method.

3. An Overlapping Schwarz Method

We describe the method in a simple case of two subdomains for a Poisson equation with a Dirichlet boundary condition as follows:∇²𝑢(𝑢=)=𝑓,inΩ,𝑢,onΓ.(3.1) Suppose Ω is split into two overlapping subregions Ω₁ and Ω₂ enclosed by 𝜕Ω₁ and 𝜕Ω₂, respectively, where 𝜕Ω_𝑖=𝜕Ω_𝑖⧵Γ_𝑖∪Γ_𝑖, 𝑖=1,2, 𝜕Ω_𝑖⧵Γ_𝑖 is the contribution of the real boundary to subdomain 𝑖 and Γ_𝑖 is the interior boundary associated with Ω_𝑖 artificially defined (Figure 2). To explain the method, let 𝑢^𝑛𝑖 denote the approximate solution on Ω_𝑖 after 𝑛 iterations, 𝑢^𝑛1|_Γ2 the restriction of 𝑢^𝑛1 to Γ₂, and 𝑢^𝑛2|_Γ1 the restriction of 𝑢^𝑛2 to Γ₁. To start the iteration, an initial approximation 𝑢⁰², for the solution in Ω₂, is required. The Schwarz method is then proceeded with the following iterations until a solution with desired accuracy is achieved.(1)Solve (3.1) for the subdomain Ω₁ as follows: 𝑢_𝑛1∇²𝑢=𝑓inΩ₁,𝑢_𝑛1=𝑢on𝜕Ω₁⧵Γ₁,=𝑢^2𝑛−1||_Γ1onΓ₁.(3.2)(2)Solve (3.1) for the subdomain Ω₂ as follows: 𝑢_𝑛2∇²𝑢=𝑓inΩ₂,𝑢_𝑛2=𝑢on𝜕Ω₂⧵Γ₂,=𝑢^𝑛1||_Γ2onΓ₂.(3.3)

Figure 2 

The domain Ω has been split into two overlapping subdomains Ω1 and Ω2 with interior boundaries Γ1 and Γ2.

(a)

(b)

(c)

(d)

Figure 3 

The adaptive partitioning algorithm is illustrated in 4 stages.

It should be noted that the values 𝑢^𝑛2|_Γ1 and 𝑢^𝑛1|_Γ2 are approximated by interpolating in Ω₂ and Ω₁, respectively. It is easy to show that performing the above iterations is equivalent to the block Gauss-Seidle method (see [20]).

4. Adapting Domain Partitioning

For the ease of description, we present the method for the case of two subdomains using a solution function 𝑢(𝑥,𝑦)=𝑒^4𝑥+𝑒^4𝑦 in the square Ω={(𝑥,𝑦)∣0≤𝑥≤1,0≤𝑦≤1} which is like the boundary layer problems. We start with a uniform mesh in Ω={(𝑥,𝑦)∣𝑎≤𝑥≤𝑏,𝑐≤𝑦≤𝑑}(4.1) (see Figure 3(a)).

Suppose Ω is split into two overlapping subdomains Ω₁ and Ω₂ with the same area. The artificial boundaries are selected such that their nodes coincide with the initial uniform mesh points used. The mesh points are equally allocated to each subregion (see Figure 3(b)). The above uniform mesh points are first distributed adaptively through the whole domain using an appropriate adaptive method with a suitable monitor function. We use the method described in Section 2 with the arc-length monitor which is appropriate for rectangular regions (see [20]).

Having applied the adaptive mesh technique, the initial equally spaced nodes will move through the domain (see Figure 3(c)). In particular, the artificial boundary nodes, which are interior to the domain, move around and take new positions which determine a new partitioning for Ω (see Figure 3(d)). More precisely, since, according to property 1, the adaptive method keeps the number of the initial mesh points, the new position of the artificial boundary nodes determines subregions Ω^1 and Ω^2 with the same number of nodes. We show this for the left subdomain. Let (𝑥₁,𝑦₁),(𝑥₂,𝑦₂),…,(𝑥_𝑛,𝑦_𝑛) be the interior boundary points and (𝑥_𝑖𝑗,𝑦_𝑖𝑗),𝑗=1,…,𝑚,𝑖=1,…,𝑛 be the nodes interior to Ω₁ for which 𝑎<𝑥_𝑖𝑗<𝑥_𝑖 and 𝑦_𝑖𝑗=𝑦_𝑖,𝑐<𝑦_𝑖<𝑑. In addition, let (𝑥^𝑖,𝑦^𝑗) and (𝑥^𝑖𝑗,𝑦^𝑖𝑗) be the adaptive points corresponding to (𝑥_𝑖,𝑦_𝑖) and (𝑥_𝑖𝑗,𝑦_𝑖𝑗), respectively, then according to property 2, we have 𝑎<𝑥^𝑖𝑗<𝑥^𝑖,𝑗=1,…,𝑚,𝑖=1,…,𝑛.

It should be noted that this partitioning, itself, roughly equidistributes the monitor function, taking into account that, on the one hand, the adaptive mesh produced in Figure 3(c) satisfies the equidistribution condition and, on the other hand, the whole mesh points are equally divided between the subdomains.

A similar algorithm is suggested for the case of more subdomains. More precisely, having produced the adaptive nodes in the single domain Ω, they are partitioned based on a partitioning of the single domain into an arbitrary number of subregions with the same area and the same number of nodes. For example, in Figure 4, an adaptive partitioning into 4 subdomains for a solution function 𝑢=𝑒^𝑐𝑥+𝑒^𝑐𝑦, is displayed for different values of 𝑐 demonstrating different amount of variations.

(a)

(b)

(c)

(d)

Figure 4 

The adaptive partitioning in the case of 4 subdomains produced for the function 𝑢=𝑒𝑐𝑥+𝑒𝑐𝑦 are displayed for various values of 𝑐.

It is worth noting that the adaptive mesh points could be applied to each subproblem separately. But this would result in a node distribution which could be only locally adapted. Since the nodes could only move through their own subdomains, the number of nodes allocating to some subregions may not be sufficient.

5. Collocation Meshless Method

In this section, we first introduce the RBFs and then describe their application to the numerical solution of PDEs based on the collocation method.

5.1. Radial Basis Functions

RBFs are known as the natural extensions of splines to multivariate interpolation. Suppose the set of points𝑥_𝑖∈Ω∣𝑖=1,2,…,𝑁(5.1) is given, where Ω is a bounded domain in 𝑅^𝑛. The radial function 𝜑∶Ω→𝑅 is used to construct the approximate function𝑠(𝑥)=^𝑁_𝑘=1𝛼_𝑘𝜑‖‖𝑥−𝑥_𝑘‖‖,(5.2) which interpolates an unknown function 𝑓 whose values at {𝑥_𝑖}^𝑁𝑖=1 are known. ‖⋅‖ represents the Euclidean norm. The unknown coefficients 𝛼_𝑘 are determined such that the following 𝑁 interpolation conditions are satisfied:𝑓𝑥_𝑖=𝑠𝑥_𝑖=^𝑁_𝑘=1𝛼_𝑘𝜑‖‖𝑥_𝑖−𝑥_𝑘‖‖,𝑖=1,…,𝑁.(5.3) There is a large class of interpolating RBFs [21] that can be used in meshless methods. These include the linear 1+𝑟, the polynomial 𝑃_𝑘(𝑟), the thin plate spline (TPS) 𝑟²log𝑟, the Gaussian exp(−𝑟²/𝛽²), and the multiquadrics √𝛽²+𝑟² (with 𝛽 a constant parameter). In this paper, we employ a generalized TPS, that is, 𝑟⁴log𝑟 which is a particular case of 𝑟^2𝑘log𝑟 (𝑘=2). These RBFs, augmented by some polynomials, are known as the natural extension of the cubic splines to the case of 2D. In fact, they are obtained by minimizing an 𝐻^𝑚 seminorm over all interpolants for which the seminorm exists. The theoretical discussion of these RBFs has been presented in [22].

5.2. Collocation Meshless Method

We describe the collocation method for a general case of PDEs in the form𝐿𝑢=𝐹,(5.4) where 𝐿=[𝐿₁,…,𝐿_𝑁]^𝑇 represents a vector of linear operations and 𝐹=[𝑓₁,…,𝑓_𝑁]^𝑇 denotes a vector containing the right-hand sides of the equations. For instance, Poisson’s equation with a Dirichlet boundary condition𝑢Δ=𝑢=𝑔,𝑓,onin𝜕ΩΩ,,(5.5) is a very simple case of equation (5.4), where 𝐿=[Δ,𝐼]^𝑇, 𝐹=[𝑓,𝑔]^𝑇 and the operators Δ and 𝐼 act on the domain Ω and the boundary 𝜕Ω, respectively.

The collocation method is used simply to express the unknown function 𝑢 in terms of the RBFs as𝑢(𝑥)=^𝑁_𝑘=1𝛼_𝑘𝜑‖‖𝑥−𝑥_𝑘‖‖(5.6) and determine the unknowns 𝛼_𝑘 in such a way that (5.6) satisfies equation (5.4) for all interpolation points. Substituting (5.6) in equation (5.4) and imposing the 𝑁 essential conditions of the collocation method lead to a linear system of equations whose coefficient matrix consists of 𝑁 row blocks, the entries of which are of the form𝐴^𝜇𝑖𝑗=𝐿_𝜇𝜑‖‖𝑥−𝑥_𝑗‖‖||_{𝑥=𝑥𝑖},𝑖=1,…,𝑁_𝜇,𝑗=1,…,𝑁,(5.7) where 𝑁_𝜇 indicates the number of nodes associated with the operator 𝐿_𝜇.

The above collocation method is referred to as a nonsymmetric collocation method due to the nonsymmetric coefficient matrices. The invertibility of the coefficient matrix can not be guaranteed [23], although in most cases a nonsingular matrix is expected. To tackle this difficulty, the symmetric collocation method, motivated by Hermitian interpolation, has been suggested. This method leads to symmetric matrices and proves the nonsingularity of interpolation matrices at the expense of double acting the operators on the RBFs [23]. Since this work is not concerned with the singularity of the matrices, the non-symmetric case will be implemented. Of course, the main technique discussed in this paper can be applied to the other case as well.

6. Numerical Results

In this section, we solve some PDEs by applying the above collocation method using a single domain and the adaptive multidomain decomposition. The results are compared with each other to show the effectiveness of the new proposed method. In each case, the method is used with some equally spaced nodes and adaptive nodes generated with various number of subdomains. The results are compared in four cases: (i) single domain with uniform mesh points, (ii) single domain with adaptive mesh points, (iii) uniform multidomain, and (iv) adaptive multidomain in the cases of 2 and 4 subdomains. A generalized TPS, 𝜙(𝑟)=𝑟⁴log𝑟, is employed as a basis function. In each example, a root mean square (RMS) error at 𝑚 collocation points is evaluated by RMSerror=⎷1𝑚^𝑚_𝑖=1𝑢_apr,𝑖−𝑢_ex,𝑖²,(6.1) where 𝑢_apr,𝑖 and 𝑢_ex,𝑖 denote the approximate and exact values of 𝑢, respectively, at a point 𝑖.

Also the arc-length is used, as a monitor function, for the mesh generation technique. In order to show the good performance of the proposed method, we consider some PDEs whose solutions are similar to those of the boundary layer problems.

Example 6.1. We apply the new method to the following PDE: Δ𝑢𝑢(𝑥=,𝑦𝑐)²=(𝑒𝑒^−𝑐𝑥_−𝑐𝑥+𝑒^−𝑐𝑦),inΩ,+𝑒^−𝑐𝑦,on𝜕Ω,(6.2) in the square Ω={(𝑥,𝑦)∣0≤𝑥,𝑦≤1} in the cases: 𝑐=20 and 𝑐=40. The exact solution for (6.2) is given by 𝑢(𝑥,𝑦)=𝑒^−𝑐𝑥+𝑒^−𝑐𝑦.

The plot of the solution function with the illustration of adaptive nodes in a single domain and multidomain in the case of 2 and 4 subdomains are displayed in Figure 5.

(a)

(b)

(c)

(d)

Figure 5 

A single domain with adaptive nodes and adaptive partitioning in the cases of 2 and 4 subdomains for the solution of Example 6.1.

The numerical errors produced by the RBF collocation method are given in Table 1 for different cases. The following notations are used in Tables 1 and 2. U. S. D:    Uniform single domain A. S. D:    Adaptive single domain U. 2 Sub.:   Uniform 2 subdomains A. 2 Sub.:   Adaptive 2 subdomains U. 4 Sub.:   Uniform 4 subdomains A. 4 Sub.:   Adaptive 4 subdomains

Table 1 

Error values for Example 6.1.

Table 2 

Error values for Example 6.2.

In the next example, we apply the new method to a nonlinear problem.

Example 6.2. One has𝑢Δ(𝑥𝑢,𝑦=)𝑓=(𝑔𝑥,(𝑥𝑦,,𝑦𝑢)),,ionnΩ𝜕,Ω,(6.3) with Ω={(𝑥,𝑦)∣0.01≤𝑥,𝑦≤1}, in two cases.
Case 1. 𝑓(𝑥,𝑦,𝑢)=4𝑢³ and 𝑔(𝑥,𝑦)=1/(𝑥+𝑦) with the exact solution 𝑢=1/(𝑥+𝑦).
Case 2. 𝑓(𝑥,𝑦,𝑢)=12𝑢² and 𝑔(𝑥,𝑦)=1/(𝑥+𝑦)² with the exact solution 𝑢=1/(𝑥+𝑦)².

The plot of the solution function with the illustration of adaptive nodes in a single domain and multidomain in the cases of 2 and 4 subdomains are displayed in Figure 6. The numerical errors produced by the RBF collocation method are given in Table 2 for different cases.

(a)

(b)

(c)

(d)

Figure 6 

A single domain with adaptive nodes and adaptive partitioning in the cases of 2 and 4 subdomains for the solution of Example 6.2.

We remark that a primary solution was first obtained by a uniform single domain. Secondly, this solution was used in approximating the derivative of the solution which was required to evaluate the monitor function in (2.2) at each collocation node. Thirdly, the primary solution was used as the initial approximate values of the interior boundary points required in the iterative Schwarz method for solving uniform DDM and adaptive DDM. Using these relatively accurate initial solutions considerably reduces the number of iterations in the Schwarz method (2-3 iterations).

As is observed, in both examples, the distribution of the mesh points is densely close to the origin due to the large variation of the solutions, in this area. This results in a small subdomain around the origin and larger ones farther from the origin. In all cases, except the latter case of Example 6.2, the numerical errors demonstrate improvements in the case of uniform DDM, a single domain with adaptive nodes and adaptive DDM. However, the best results are achieved in the case of using the new method, adaptive domain decomposition method. This is due to the fact that, although using a single domain with adaptive nodes improve the accuracy of the solution, the interpolation matrix, in this case, may be a large matrix with a large condition number. Consequently, applying the adaptive DDM improves the conditioning while taking advantage of the equidistribution property which leads to more accurate solution.

In the latter case of Example 6.2, the variation of the solution is much larger than that of the former case. As a consequence, the distribution of the nodes is extensively dense around the origin. This, itself, causes the conditioning of the problem and leads to less accurate solution. But, the numerical results are improved in the case of applying the adaptive DDM. This improvement is justified by the fact that, in the case of the adaptive DDM, the interpolation matrix is split into some smaller matrices which improve the condition numbers of the local matrices.The results show that application of the uniform DDM may not be adequate for these sort of solution functions.

In addition to improvement of the accuracy, the computational efficiency is considerably improved due to the use of the small size matrices in the subproblems.

7. Conclusion

A new domain partitioning technique based on adaptive nodes was presented to improve meshless type methods. The method was based on first, generating some adaptive nodes in a rectangle and then partitioning the domain based on the distribution of the points.

The method was applied to the collocation meshless method; however, it is applicable to any sort of mesh-free method in which the connectivity of the nodes is not important. The method was applied to the case of rectangular domains using an adaptive technique based on equi-distributing principle. One may be interested in adapting the proposed method with the other mesh generation technique. Applying the new method to the irregular regions is more complicated and is currently being considered. The new method was examined by considering some PDEs solved by a collocation meshless method and the results demonstrated considerable reduction in the error values.