Enriched Meshfree Method for an Accurate Numerical Solution of the Motz Problem
We present an enriched meshfree solution of the Motz problem. The Motz problem has been known as a benchmark problem to verify the efficiency of numerical methods in the presence of a jump boundary data singularity at a point, where an abrupt change occurs for the boundary condition. We propose a singular basis function enrichment technique in the context of partition of unity based meshfree method. We take the leading terms of the local series expansion at the point singularity and use them as enrichment functions for the local approximation space. As a result, we obtain highly accurate leading coefficients of the Motz problem that are comparable to the most accurate numerical solution. The proposed singular enrichment technique is highly effective in the case of the local series expansion of the solution being known. The enrichment technique that is used in this study can be applied to monotone singularities (of type with ) as well as oscillating singularities (of type ). It is the first attempt to apply singular meshfree enrichment technique to the Motz problem.
The Motz problem was first introduced by Motz . Since it was introduced, the Motz problem has served as a benchmark problem to verify the efficiency of numerical methods in the presence of a singularity. The Motz problem is a Laplace equation with mixed Neumann-Dirichlet boundary conditions as follows:
The computational domain is shown in Figure 1. The Motz problem has a jump boundary data singularity at the origin where the boundary condition abruptly changes from the Neumann boundary condition to the Dirichlet boundary condition. The exact solution of the Motz problem can be obtained by conformal mapping .
On the other hand, the asymptotic solution of (1) can be expressed as an infinite series,and the series expansion is valid on the entire domain [3, 4]. Therefore, the following partial sum of (2) is called an approximate solution of the Motz problem:
Remark 1. The function for each satisfies the Dirichlet boundary condition along the negative -axis and the Neumann boundary condition along the positive -axis.
There have been many efforts to get an accurate estimate for the singular coefficients in (3). Special finite difference method  and global element methods [6, 7] were introduced to determine . Other methods such as one-zone blending , combination of two-zone blended singular basis functions , boundary method , the method of auxiliary mapping (MAM), and mesh refinement in the vicinity of singularity  are available for this problem. However, approximating is generally unsatisfactory or hard to compute. On the other hand, a finite difference approach , collocation Trefftz method , and reproducing polynomial boundary particle methods (RPBPM)  successfully estimated the first several coefficients. The collocation Trefftz method and RPBPM especially give a highly accurate leading coefficients for the Motz problem.
In recent years many meshless methods showed great success: Element-Free Galerkin Method (EFGM) , Method of Finite Sphere , Reproducing Kernel Particle Method (RKPM) , and Reproducing Polynomial Particle Method (RPPM) [17, 18] to name but a few. These methods do not explicitly use a traditional mesh. On the other hand, there is partition of unity methods that use finite element mesh explicitly. Thus, in essence it is not purely meshfree. However, it is widely known as a meshfree method because the usage of a background mesh is completely different from the conventional finite element method. For these methods, the mesh is used to construct a partition of unity rather than to build shape functions as in the conventional finite element method. h-p clouds , Extended Finite Element Method (XFEM) [20, 21], Generalized Finite Element Method (GFEM) [22–24], and mesh based construction of flat-top partition of unity functions (MFPUM)  are in this category. Although the usage of a mesh is completely different for meshfree method and conventional FEM, there has been an effort to couple these methods on a single conventional mesh .
Using only the piecewise polynomial partition of unity shape functions as a basis function for the meshfree approximation can not reproduce polynomials other than a constant. Therefore, it is important to choose an appropriate local approximation function that has higher polynomial reproducing order. In addition, if the local approximation functions have the Kronecker delta property, imposing essential boundary conditions becomes simple as conventional FEM. On top of that one can create highly regular basis functions that have compact support. This enables us to solve higher order partial differential equations such as biharmonic or polyharmonic problems [27–29] without using the Hermite finite elements, which are difficult to implement in higher dimension. A 3D extension of the partition of unity based meshfree method, which uses compactly supported highly regular basis functions that have the Kronecker delta property, is also available .
Many meshfree methods [4, 16, 20, 21] have been successfully applied to engineering problems, especially crack modeling. An enrichment with discontinuous functions was the main interest since the discontinuity naturally arises in the solution of the elasticity problem with a cracked domain. Likewise, it is natural to enrich basis functions that can closely approximate the solution of the Motz problem. We use meshfree method with singular enrichment functions, which are taken from the series expansion. The piecewise polynomial partition of unity function helps us to enrich the local approximation space with singular functions. As a consequence, we obtain highly accurate approximation for the singular coefficients as part of our meshfree solution. The enrichment technique that is introduced in this study can be used to solve problems that have a known local series expansion or even an asymptotic series expansion.
2. An Enriched Meshfree Method
The corresponding variational equation for the Motz problem, given in (1), is as follows.
Find such that on andfor all
Now, the meshfree method based on the Generalized Finite Element Method for a numerical solution to this problem is described as follows.(1)Generate a Background Mesh. To construct PU functions with flat-top and highly smooth local approximation functions for numerical solutions of (4), the domain is partitioned into polygonal patches. We use quadrangular patches for convenience. Unlike the conventional FEM mesh, we allow the hanging nodes for the background mesh.(2)Construct Partition of Unity Functions with a Flat-Top. For , let be the PU function with flat-top corresponding to the patches and let be the support of .(3)Planting Particles. Let be the reference rectangular patch. Suppose , are arbitrarily (or uniformly) distributed particles on the reference patches. Also, let (if the physical patch is quadrangular) be the patch mappings. Then, through the patch mappings , we have particles , on .(4)Local Approximation Functions. Suppose , are smooth shape functions corresponding to the particles planted in the reference patch that satisfy the Kronecker delta property. Then these shape functions on the reference patch can be used to build local approximation functions on the support of the physical patch as follows: Here, we use the tensor product of Lagrange interpolation functions corresponding to arbitrarily spaced numbers of nodes in the -direction and those corresponding to arbitrarily spaced numbers of nodes in the -direction. These shape functions satisfy the Kronecker delta property. In other words, let be the th Lagrange interpolation polynomial of degree associated with distinct nodes in , defined by Then are shape functions with polynomial reproducing order .(5)Smooth Global Basis Functions. The global approximation functions with compact support are constructed as follows: These global approximation functions are highly smooth and correspond to the particles: We also assume that each set of local approximation functions has the polynomial reproducing property and satisfies the Kronecker delta property at the particles planted on .(6)Meshfree Approximation Space. The vector space spanned by those approximation functions defined by (7), denoted by , is called the meshfree approximation space. Then the meshfree approximation is given by The meshfree approximation is essentially a Galerkin approximation with the use of . Find such that
2.1. Basis Enrichment
There are many different choices for the local approximation functions. Some of them can be found in . It is important to understand that this flexibility of selecting local approximation functions is one of the most powerful aspects of partition of unity based meshfree methods. The optimal choice of local approximation spaces is discussed extensively in .
For the Motz problem, we use singular functions on the patch, which includes the point singularity, and on the rest of the patches we use Lagrange interpolation polynomials for the local approximation functions. In case of using singular basis functions, the numerical integration of the bilinear form will not be exact. As pointed out earlier, the flexibility to choose a local approximation space enables us to add special functions to approximation spaces. Enrichment means adding some special functions to the sets of existing local approximation functions in the following way.
Definition 2. The enriched meshfree approximation space is defined as follows. Suppose one wants to enrich local approximation functions on the patches with sets of special functions: Then is the vector space spanned by
The enriched meshfree method for the second-order elliptic equation with the boundary condition (1) associated with the enriched meshfree approximation space is as follows.
Findsuch thatwhere and are defined as (4).
Calculation of a stiffness matrix involves the following integration:where and indicate the patch number and . Also and the components of the load vector are given as follows:
Remark 3. Since the PU function is piecewise polynomial, if the local approximation functions are polynomials as well, then integration (16) can be calculated exactly. Note that integrations, (16) and (17), contain nonpolynomial singular function. In this case we need to approximate the numerical integration.
Integration (16) has four different types:where is a regular piecewise polynomial approximation function and is a singular enrichment function.
2.2. Error Estimate
The error estimate of partition of unity based methods can be found in several different forms. We give error estimate based on .
Theorem 4 (let be a 2D convex polygon). Also let be a partition of unity corresponding to the patch of . Also denote for the flat-top region of . Assume there exists a positive integer M such that, for every , , where . Let a collection of local approximation spaces be given on each patch .
Let us assume the local approximation space has the following approximation properties on each patch .
The function can be approximated by a function such that (1),(2),(3),for all Then, satisfies the following global estimates:where is a constant independent of .
The most popular method to solve singularity problems by the conventional FEM is using adaptive mesh refinement. Since hanging nodes are not allowed in the traditional FEM, the mesh has to be globally refined to avoid hanging nodes. Moreover, the convergence could be slower than expected even if the mesh becomes extremely fine . Our goal is to numerically calculate the system of equation that is given in (15). Since the partition of unity allows us to use hanging nodes in the background mesh, we partition into 39 patches as shown in Figure 2: 38 rectangular patches along the boundary and one large rectangular patch containing the singularity. The enriched patch is shaded in Figure 3. The large patch, which contains the point singularity, is called the enriched patch, where singular functions are used as local approximation functions. On the remaining 38 small patches, the local approximation functions defined by (7) are used. Two neighboring patches have overlapping supports for the partition of unity functions. The overlapping parts are strips of -width. In our numerical test, we use . The thin strips represent the overlapping regions in Figure 2. Except on the overlapping regions, partition of unity functions has the value of one to ensure the linear independence of the local approximation functions.
Special attention is needed for numerical integration given in (19)–(21) on the enriched patch. Since all local approximation functions on the enriched patch are singular functions, we do not have integration (18) on the enriched patch. Likewise, on other patches we do not have integration of type (21). Integrals in (19) and (20) seem unbounded because these have singular integrands. However, these integrals do exist because the integration is only performed on overlapping regions with strips along the boundary of the enriched patch. Although the numerical integration may not be exact, the integrands are smooth since the overlapping region is far away from the singular point, . We should get accurate result for these integrations with a reasonable number of quadrature points. Now numerical integration of (21) is the only integral that is left to consider. The following theorem not only shows the boundness of the integral given in (21) but also gives an idea of how to get an accurate numerical integration for the proposed method.
Theorem 5. Considerwhere denotes the flat-top region of the PU function corresponding to the enriched patch shown in Figure 3 and .
Proof. By a coordinate transformation, the gradient operator is transformed as follows: where . For simplicity, we suppress in . Then,Since , Plugging (27) into (26), we have the following: which proves the theorem.
Let be the global indices corresponding to the local enrichment functions , respectively, and let be component of stiffness matrix . Then is given by the following integration:where and are the supports of the partition of unity functions and , respectively. Recall that the functions and are globally defined functions. Then, the overlapping region can be decomposed as : a disjoint union of and . Since on , of (29) can be written: where and .
Let us examine first. Note that is a piecewise polynomial and the integral region is far away from the singularity of . Also it is important to note that the integral region has a rectangular shape. Therefore, integral is straight forward to evaluate. With enough Gauss points, integral should be nearly exact.
Let be the coordinate transformation, and, and let be the patch mapping from the reference patch to a quadrangle patch . See Figure 4(a). Then we haveSimilarly, Note , , , and are all function of . Let us define . Then, where is the number of Legendre-Gauss quadrature points, is the Legendre-Gauss point, and is the th quadrature weight.
(a) A rectangular patch mapping
(b) A blending mapping
On the other hand, by Theorem 5, Part is as follows: of the enriched patch can be further decomposed into two separate regions and . in Figure 3 is the half disk and is its complement; that is, (see Figure 3). On , the integral of part can be obtained exactly.
Then, integral over can be simplified as follows:Using (35), we have where is the Kronecker delta. Now over region, we need to perform numerical integration because of its irregular shape. Since region is arch-shaped, we need to use the patch mapping of the blending type for this integration.
Remark 6. Note that the integrand becomes highly oscillatory as the value of increases. For instance, for the integrand oscillates about 20 times when the angle moves from 0 to . Therefore, we need to divide the integral domain in -direction. For an accurate numerical integration, -direction integral was subdivided into 32 integrals to achieve an accurate numerical integration.
The blending mapping , Figure 4(b), maps the reference patch onto a quadrangle with a curved side. Let and in terms of polar coordinates and and in rectangular coordinates. The blending patch mapping is defined as follows: Also its partial derivatives are as follows: The determinant of Jacobian (), denoted by , is . Hence integral over region can be transformed over region : As pointed out earlier to obtain accurate numerical integration, region should be partitioned into many pieces to capture the oscillatory behavior of the integrand. Let us say is partitioned into pieces. Also let be the transformation from the reference square into the integral region , where . Then integral becomes the following: where denotes the number of Legendre-Gauss quadrature points, denotes the Legendre-Gauss point, and is the th quadrature weight. Thus, we have completed the calculation of . On the other hand, since the Motz problem is a Laplace equation, calculating the load vector for (15) is straightforward.
Since we have the Kronecker delta property of the basis functions, imposing essential boundary condition is as simple as in the conventional FEM. Therefore, the solution of the Motz problem does not have any error on the Dirichlet boundary conditions because the essential boundary conditions are constants and local approximation functions have the polynomial reproducing property.
3. Numerical Results
The numerical tests are carried out by using the local approximation functions that have polynomial reproducing orders 2, 4, 6, and 8 on each of the 38 rectangular patches that do not contain the singularity and 40 singular functions on the enrichment patch containing the singularity.
Table 1 contains the numerical results obtained by the proposed method, which uses singular functions enrichment for the Motz problem. It is known that the numerical solution of the Motz problem obtained by  is the most accurate one. In Table 1, the errors in the max norm as well as the energy norm are computed using the best known solution  as the true solution. The strain energy of is . The relative error in the energy norm is defined by where and are the strain energy obtained by the solution of proposed method and provided in the literature , respectively.
Figure 5 shows the error of the numerical solution. On the enrichment patch, we see that there is virtually no error near the point singularity. The singular behavior of the Motz problem was effectively captured by 40 enriched singular functions on the enriched patch and local approximation functions with polynomial reproducing order 8 on the surrounding 38 patches. Also, the error plot shows that the meshfree solution does not have any error along the Dirichlet boundary, and . It is worthwhile to note that the degrees of freedom to capture the singularity at the enrichment patch are only 40, which is the maximum number of enrichment functions we have used. We have chosen 40 enrichment functions because adding more enrichment functions does not enhance the accuracy of the numerical solution in double precision, as the coefficients of series expansion decay exponentially to zero.
If the structure of the singularity is known, RSPM (Reproducing Singularity Particle Methods) is known to be more effective than the adaptive RPPM . RSPM is a Galerkin approximation associated with the use of RSP shape functions on the patches containing singularities and with the use of RPP shape functions on other patches for local approximation functions. Thus, RSPM is similar to the method of auxiliary mapping [10, 34–36] in the framework of the -version FEM. Table 2 compares the result of the proposed method with RSPM and -FEM with MAM. We conclude that the enriched meshfree method is superior to both RSPM and -FEM with MAM for the Motz problem. Also it is important to note that both -version FEM with MAM and RSPM can not estimate in (3) but the proposed method can.
We emphasize the fact that our enriched solution has no errors along and because the local approximation functions either satisfy the boundary condition or have the polynomial reproducing and Kronecker delta property. Our numerical solution satisfies the boundary condition exactly. Hence, along , the enriched meshfree solution is clearly better than the best known solution , which is oscillatory on
Remark 7. Those existing numerical methods, listed in Table 4, such as boundary methods , collocation Trefftz methods , and Integrated Singular Basis Function Methods , need special treatment to accurately impose Dirichlet boundary condition. However, the proposed method can impose the essential boundary conditions straight forwardly using the Kronecker delta property of the local approximation functions.
Unlike Li’s solution, which has the power series form given in (3) on the entire domain, the numerical solution of the proposed method has the representation given in (14). Thus in (14) are not intended to be used in (3). However, ’s in our proposed method can be regarded as the coefficients of the power series, since we only have singular functions in the enrichment patch and the partition of unity function , in (14), has the value of one if we restrict the numerical solution locally near the singularity. Therefore, we may compare in (14) with Li’s solution , which is the best known numerical solution for the Motz problem. We denote the coefficients from our method to and Li’s solution to . , , is listed in Table 3. This time we compare the errors of and along the boundary where the true solution is 500. As mentioned earlier, the enriched meshfree solution given in (14) does not have any error on the boundary . To compare the accuracy with the best known solution, , , are extracted from the meshfree solution and applied to (3). Therefore, the series solution obtained by using , , as coefficients exhibits error on and they are measured by and , respectively. The errors along the boundary are plotted in Figure 6.
Remark 8. The absolute error in maximum norm on is measured as for the series solution with and for the series solution with . In the mean square norm, the errors are and , respectively. Although Li’s coefficients result in series solution that is slightly accurate on in both maximum and mean square norm, the proposed method has more significant digits for leading coefficients; see Table 4.
Finally to show the effectiveness of the proposed method, the first 25 coefficients of the partial sum (3) obtained by various numerical methods are compared in Table 4. We conclude that the proposed method results in highly accurate coefficients . The first coefficient especially has at least 13 significant digits.
4. Concluding Remark
We introduced how to choose and enrich singular functions in the case of the local series expansion of the solution being known for the partition of unity based meshfree method. The enrichment technique that is used in this study can be applied to not only monotone singularities (of type with ) but also oscillating singularities (of type ).
The author declares that there are no competing interests regarding the publication of this paper.
This work is supported by Gachon University, Grant no. GCU-2014-0130.
M. Motz, “The treatment of singularities of partial differential equations by relaxation methods,” Quarterly of Applied Mathematics, vol. 4, pp. 371–377, 1947.View at: Google Scholar
J. B. Rosser and N. Papamichael, “A power series solution of a harmonic mixed boundary value problem,” Technical Summary Report, 1975.View at: Google Scholar
H.-S. Oh and J. W. Jeong, “Reproducing polynomial (singularity) particle methods and adaptive meshless methods for two-dimensional elliptic boundary value problems,” Computer Methods in Applied Mechanics and Engineering, vol. 198, no. 9-12, pp. 933–946, 2009.View at: Publisher Site | Google Scholar | Zentralblatt MATH