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 basedmeshfreemethod.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 rα with α < 1) as well as oscillating singularities (of type rα sin(ε log r)). It is the first attempt to apply singular meshfree enrichment technique to the Motz problem.


Introduction
The Motz problem was first introduced by Motz [1].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 (0, 0) 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 [2].
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  +1/2 cos( + 1/2) 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 d in (3).Special finite difference method [5] and global element methods [6,7] were introduced to determine d .Other methods such as onezone blending [3], combination of two-zone blended singular basis functions [8], boundary method [9], the method of auxiliary mapping (MAM), and mesh refinement in the vicinity of singularity [10] are available for this problem.However, approximating d is generally unsatisfactory or hard to compute.On the other hand, a finite difference approach [11], collocation Trefftz method [12], and reproducing polynomial boundary particle methods (RPBPM) [13] 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) [14], Method of Finite Sphere [15], Reproducing Kernel Particle Method (RKPM) [16], 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 [19], Extended Finite Element Method (XFEM) [20,21], Generalized Finite Element Method (GFEM) [22][23][24], and mesh based construction of flat-top partition of unity functions (MFPUM) [25] 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 [26].
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][28][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 [30].
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 d 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.

An Enriched Meshfree Method
The corresponding variational equation for the Motz problem, given in (1), is as follows.
Find  ∈  1 (Ω) such that  =   on Γ  and for all V ∈  1  (Ω).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.
These global approximation functions are highly smooth and correspond to the particles: We also assume that each set {  :  = 1, . . .,   } 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  mfree , 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  mfree .
Find  * ∈  mfree such that

Basis Enrichment.
There are many different choices for the local approximation functions.Some of them can be found in [27].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 [32].
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.
Then  enrich 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  enrich is as follows. Find such that where B and F are defined as (4).Calculation of a stiffness matrix involves the following integration: where  and  indicate the patch number and ∇ = (/, /)  .Also V , , V , ∈  enrich 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.

Error Estimate.
The error estimate of partition of unity based methods can be found in several different forms.We give error estimate based on [33].
Let us assume the local approximation space   has the following approximation properties on each patch   .
The function  ∈  1 (Ω) can be approximated by a function , for all  = 1, 2, . . ., .Then, satisfies the following global estimates: where  ∇ is a constant independent of .

Implementation.
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 [10].
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 2-width.In our numerical test, we use  = 0.01.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 2 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, (0, 0).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.
where Ω flat 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  = ( cos  −(1/) sin  sin  (1/) cos  ).For simplicity, we suppress (, ) in (, ).Then, Since   (, ) =  (1/2+) cos((1/2 + )), 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 (, )th component of stiffness matrix .Then   is given by the following integration:
Let us examine () first.Note that Ψ  is a piecewise polynomial and the integral region  nonflat 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,  =  cos  and  =  sin , and let   : (, ) → (, ) be the patch mapping from the reference patch P to a quadrangle patch R. See Figure 4(a).Then we have Similarly, where   is the number of Legendre-Gauss quadrature points,   is the Legendre-Gauss point, and   is the th quadrature weight.
On the other hand, by Theorem 5, Part () is as follows: flat of the enriched patch can be further decomposed into two separate regions  and   . = {(, ) | 0 ≤  ≤ , 0 ≤  ≤ } in Figure 3 is the half disk and   is its complement; that is,  flat =  ∪   (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 cos(( − )) becomes highly oscillatory as the value of | − | increases.For instance, for  = 40,  = 1 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 P onto a quadrangle Q with a curved side.Let  1 = (,  2 ) and  2 = (,  1 ) in terms of polar coordinates and  3 = ( 3 ,  3 ) and  4 = ( 4 ,  4 ) in rectangular coordinates.The blending patch mapping   is defined as follows: )) . ( Also its partial derivatives are as follows: The 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   = ⋃  =1    .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.

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 [9] 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 [9] as the true solution.The strain energy of  ∈  1 (Ω) is U = (1/2)B(, ).The relative error in the energy norm is defined by where U * and U Li are the strain energy obtained by the solution of proposed method and provided in the literature [9], 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, Γ 1 and Γ 2 .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 [33].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][35][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.We emphasize the fact that our enriched solution has no errors along Γ 1 and Γ 2 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 Γ 2 , the enriched meshfree solution is clearly better than the best known solution [9], which is oscillatory on Γ 2 .
Remark 7.Those existing numerical methods, listed in Table 4, such as boundary methods [9], collocation Trefftz methods [12], and Integrated Singular Basis Function Methods [4], 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, d '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 [9], which is the best known numerical solution for the Motz problem.We denote the coefficients

d †
+1/2 cos( + 1/2)|, respectively.The errors along the boundary Γ 2 are plotted in Figure 6.slightly accurate on Γ 2 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 d 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 d .The first coefficient especially has at least 13 significant digits.

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  < 1) but also oscillating singularities (of type   sin( log ) [36]).

Competing Interests
The author declares that there are no competing interests regarding the publication of this paper.

Figure 1 :
Figure 1: The domain for the Motz problem.Essential boundary conditions are imposed on Γ 1 and Γ 2 ; otherwise Neumann boundary condition is imposed.

Figure 2 :
Figure 2: 39 overlapping patches, which consist of one large patch and 38 surrounding patches, cover the domain.Each patch is the support of corresponding partition of unity function.The overlapping regions are pictured as strip that has width 2.

Figure 3 :
Figure 3: Shaded rectangle is the enriched patch, which is subdivided by disjoint regions  and   for numerical integration.

Figure 5 :
Figure 5: Maximum error plot when local approximation functions have polynomial reproducing order 8 for the 38 surrounding patches and 40 singular basis functions are used on enriched patch.

Table 2 :
Relative errors are compared for different  degrees in energy norm.Also it is important to note that both -version FEM with MAM and RSPM can not estimate d in (3) but the proposed method can.

Table 3 :
Computed coefficients d * of ∑ The absolute error in maximum norm on Γ 2 is measured as 6.89 × 10 −9 for the series solution with d *  and 5.47 × 10 −9 for the series solution with d †  .In the mean square norm, the errors are 2.78 × 10 −9 and 2.26 × 10 −9 , respectively.Although Li's coefficients result in series solution that is

Table 4 :
Comparison of the first 25 coefficients obtained by other methods.