An Adaptive Domain Partitioning Technique for Meshfree-type Methods

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.


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.

Adaptive Nodes in a Rectangle
For a 2D domain Ω, a partition {Ω i } n i 1 is called equidistributing with respect to a nonnegative piecewise continuous function M x, y see 17 if where the monitor function M 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

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 The equidistribution process is performed in three stages.In the first stage, equidistribution is performed in the direction of the x-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 x coordinate direction again Figure 1 c .The mesh resulting from the above procedure forms quadrilaterals whose sides equidistribute 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 is the initial set of equally spaced nodes and is the set of points produced by the adapting strategy, then N N and there will be a one-to-one correspondence between the members of A and B. In fact, the adaptive method moves a point x i , y i to a new position denoted by x i , y i .
2 It keeps the order of points, that is, for every two points x i , y i and x j , y j , if x i < x j , then x i < x j , and if y i < y j , then y i < y j .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.

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: u u, on Γ.

3.1
Suppose Ω is split into two overlapping subregions Ω 1 and Ω 2 enclosed by ∂Ω 1 and ∂Ω 2 , respectively, where after n iterations, u n 1 | Γ 2 the restriction of u n 1 to Γ 2 , and u n 2 | Γ 1 the restriction of u n 2 to Γ 1 .To start the iteration, an initial approximation u 0 2 , for the solution in Ω 2 , 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 Ω 1 as follows: 3.2 2 Solve 3.1 for the subdomain Ω 2 as follows:

3.3
It should be noted that the values u n 2 | Γ 1 and u n 1 | Γ 2 are approximated by interpolating in Ω 2 and Ω 1 , respectively.It is easy to show that performing the above iterations is equivalent to the block Gauss-Seidle method see 20 .

Adapting Domain Partitioning
For the ease of description, we present the method for the case of two subdomains using a solution function u x, y e 4x e 4y in the square Ω { x, y | 0 ≤ x ≤ 1, 0 ≤ y ≤ 1} which is like the boundary layer problems.We start with a uniform mesh in see Figure 3 a .Suppose Ω is split into two overlapping subdomains Ω 1 and Ω 2 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 x 1 , y 1 , x 2 , y 2 , . . ., x n , y n be the interior boundary points and x ij , y ij , j 1, . . ., m, i 1, . . ., n be the nodes interior to Ω 1 for which a < x ij < x i and y ij y i , c < y i < d.In addition, let x i , y j and x ij , y ij be the adaptive points corresponding to x i , y i and x ij , y ij , respectively, then according to property 2, we have a < x ij < x i , j 1, . . ., m, i 1, . . ., n.
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 u e cx e cy , is displayed for different values of c demonstrating different amount of variations.
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.

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.

Radial Basis Functions
RBFs are known as the natural extensions of splines to multivariate interpolation.Suppose the set of points is given, where Ω is a bounded domain in R n .The radial function ϕ : Ω → R is used to construct the approximate function which interpolates an unknown function f whose values at {x i } N i 1 are known.• represents the Euclidean norm.The unknown coefficients α k are determined such that the following N interpolation conditions are satisfied: There is a large class of interpolating RBFs 21 that can be used in meshless methods.These include the linear 1 r, the polynomial P k r , the thin plate spline TPS r 2 log r, the Gaussian exp −r 2 /β 2 , and the multiquadrics β 2 r 2 with β a constant parameter .In this paper, we employ a generalized TPS, that is, r 4 log r which is a particular case of r 2k log r k 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 H m seminorm over all interpolants for which the seminorm exists.The theoretical discussion of these RBFs has been presented in 22 .

Collocation Meshless Method
We describe the collocation method for a general case of PDEs in the form The collocation method is used simply to express the unknown function u in terms of the RBFs as and determine the unknowns α k 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 N essential conditions of the collocation method lead to a linear system of equations whose coefficient matrix consists of N row blocks, the entries of which are of the form where N μ indicates the number of nodes associated with the operator L μ .
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.

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, φ r r 4 log r, is employed as a basis function.In each example, a root mean square RMS error at m collocation points is evaluated by where u apr,i and u ex,i denote the approximate and exact values of u, respectively, at a point i.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: 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.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  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.
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.

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.

Figure 1 :
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.

Figure 3 :
Figure 3: The adaptive partitioning algorithm is illustrated in 4 stages.

Adaptive partition c = 2 b partition c = 4 cpartition c = 6 dFigure 4 :
Figure 4: The adaptive partitioning in the case of 4 subdomains produced for the function u e cx e cy are displayed for various values of c.

1
, . . ., L N T represents a vector of linear operations and F f 1 , . . ., f N T denotes a vector containing the right-hand sides of the equations.For instance, Poisson's equation with a Dirichlet boundary condition Δu f, in Ω, u g, on ∂Ω, 5.5 is a very simple case of equation 5.4 , where L Δ, I T , F f, g T and the operators Δ and I act on the domain Ω and the boundary ∂Ω, respectively.
u x, y e −cx e −cy , on ∂Ω, 6.2 in the square Ω { x, y | 0 ≤ x, y ≤ 1} in the cases: c 20 and c 40.The exact solution for 6.2 is given by u x, y e −cx e −cy .

Figure 5 :
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.

Figure 6 :
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. .

Table 1 :
Error values for Example 6.1.

Table 2 :
Error values for Example 6.2.