# Enriched Meshfree Method for an Accurate Numerical Solution of the Motz Problem

Authors: Won-Tak Hong
Journal: Advances in Mathematical Physics
 (2016)

---

## Abstract

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


---

## Body

## 1. 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:(1)
-
Δ
u
=
0
in
Ω
=
x
,
y
∣
-
1
<
x
<
1
,
0
<
y
<
1
,
u
=
0
on
Γ
1
,
u
=
500
on
Γ
2
,
∂
u
∂
n
=
0
at
Γ
∖
Γ
1
∪
Γ
2
.The computational domain is shown in Figure1. 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].Figure 1
The domain for the Motz problem. Essential boundary conditions are imposed onΓ
1 and Γ
2; otherwise Neumann boundary condition is imposed.On the other hand, the asymptotic solution of (1) can be expressed as an infinite series,(2)
u
r
,
θ
=
∑
i
=
0
∞
d
i
r
i
+
1
/
2
cos
⁡
i
+
1
2
θ
,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:(3)
u
N
r
,
θ
=
∑
i
=
0
N
d
~
i
r
i
+
1
/
2
cos
⁡
i
+
1
2
θ
.Remark 1.
The functionr
i
+
1
/
2
cos
⁡
(
i
+
1
/
2
)
θ for each i satisfies the Dirichlet boundary condition along the negative x-axis and the Neumann boundary condition along the positive x-axis.There have been many efforts to get an accurate estimate for the singular coefficientsd
~
i in (3). Special finite difference method [5] and global element methods [6, 7] were introduced to determine d
~
i. Other methods such as one-zone 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
~
i 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–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–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
~
i 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.Findu
∈
H
1
(
Ω
) such that u
=
u
d on Γ
D and(4)
B
u
,
v
≡
∫
Ω
∇
u
·
∇
v
-
∫
Γ
D
v
∇
u
·
n
=
∫
Ω
f
v
+
∫
Γ
N
u
n
v
≡
F
v
,for all v
∈
H
D
1
(
Ω
)
.Now, the meshfree method based on theGeneralized 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 J
=
1,2
,
…
,
M, let Ψ
J be the PU function with flat-top corresponding to the patches Q
J and let ω
J be the support of Ψ
J.
(3)
Planting Particles. Let Q
r be the reference rectangular patch. Suppose p
^
k
,
k
=
1,2
,
…
,
N, are arbitrarily (or uniformly) distributed particles on the reference patches. Also, let ϕ
J
q
:
Q
r
→
ω
J (if the physical patch is quadrangular) be the patch mappings. Then, through the patch mappings ϕ
J
q, we have particles p
J
k
=
(
ϕ
J
q
)
-
1
(
p
^
k
), k
=
1
,
…
,
N
J on ω
J.
(4)
Local Approximation Functions. Suppose g
^
k
,
k
=
1
,
…
,
N
^, 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 w
J of the physical patch Q
J as follows:(5)
g
J
k
=
g
^
k
∘
ϕ
J
q
-
1
,
k
=
1
,
…
,
N
J
.
Here, we use the tensor product of Lagrange interpolation functions corresponding to arbitrarily spacedN
x numbers of nodes in the x-direction and those corresponding to arbitrarily spaced N
y numbers of nodes in the y-direction. These shape functions satisfy the Kronecker delta property. In other words, let L
n
,
j
(
x
) be the jth Lagrange interpolation polynomial of degree (
n
-
1
) associated with n distinct nodes x
1
,
…
,
x
n in [
0,1
], defined by(6)
L
n
,
j
x
=
∏
k
=
1
,
k
≠
j
n
x
-
x
k
∏
k
=
1
,
k
≠
j
n
x
j
-
x
k
.
Theng
J
k
(
x
,
y
)
=
L
n
,
l
(
x
)
×
L
n
,
m
(
y
)
;
1
≤
l
,
m
≤
n
,
k
=
n
(
l
-
1
)
+
m are shape functions with polynomial reproducing order n
-
1.
(5)
Smooth Global Basis Functions. The global approximation functions with compact support are constructed as follows:(7)
Φ
J
k
x
,
y
=
Ψ
J
x
,
y
·
g
J
k
x
,
y
,
J
=
1
,
…
,
M
;
k
=
1
,
…
,
N
J
.
These global approximation functions arehighly smooth and correspond to the particles:(8)
p
J
k
,
J
=
1,2
,
…
,
M
;
k
=
1,2
,
…
,
N
J
.
We also assume that each set{
g
J
k
:
k
=
1
,
…
,
N
J
} of local approximation functions has the polynomial reproducing property and satisfies the Kronecker delta property at the particles planted on ω
J.
(6)
Meshfree Approximation Space. The vector space spanned by those approximation functions defined by (7), denoted by V
mfree, is called the meshfree approximation space. Then the meshfree approximation is given by(9)
u
∗
x
,
y
=
∑
J
∑
k
c
J
k
·
G
J
k
x
,
y
=
∑
J
Ψ
J
x
,
y
∑
k
c
J
k
·
g
J
k
x
,
y
.
The meshfree approximation is essentially a Galerkin approximation with the use ofV
mfree.
Findu
∗
∈
V
mfree such that(10)
B
u
∗
,
v
=
F
v
∀
v
∈
V
mfree
.


### 2.1. 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.Definition 2.
The enriched meshfree approximation spaceV
enrich is defined as follows. Suppose one wants to enrich local approximation functions on the patches (11)
Q
J
1
,
Q
J
2
,
…
,
Q
J
swith sets of special functions: (12)
σ
J
i
,
k
,
k
=
1,2
,
…
,
N
J
i
for
i
=
1,2
,
…
,
s
.Then V
enrich is the vector space spanned by (13)
V
mfree
∪
Ψ
J
i
σ
J
i
k
:
k
=
1,2
,
…
,
N
J
i
,
i
=
1,2
,
…
,
s
.The enriched meshfree method for the second-order elliptic equation with the boundary condition (1) associated with the enriched meshfree approximation space V
enrich is as follows.Find(14)
u
enrich
=
∑
J
=
1
M
∑
k
=
1
N
J
c
J
k
Ψ
J
g
J
k
+
∑
i
=
1
s
∑
l
=
1
N
J
i
d
J
i
,
l
Ψ
J
i
σ
J
i
,
l
,such that(15)
B
u
enrich
,
v
=
F
v
,
for
v
∈
V
enrich
,where B and F are defined as (4).Calculation of a stiffness matrix involves the following integration:(16)
∫
ω
J
∩
ω
L
∇
v
J
,
i
T
·
∇
v
L
,
j
d
Ω
,where J and L indicate the patch number and ∇
=
(
∂
/
∂
x
,
∂
/
∂
y
)
T. Also v
J
,
i
,
v
L
,
j
∈
V
enrich and the components of the load vector are given as follows:(17)
∫
ω
J
∩
ω
L
f
·
v
J
,
k
d
Ω
.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:(18)
∫
ω
J
∩
ω
L
∇
Ψ
J
·
g
J
k
T
·
∇
Ψ
L
·
g
L
l
d
Ω
,
(19)
∫
ω
J
∩
ω
L
∇
Ψ
J
·
g
J
k
T
·
∇
Ψ
L
·
σ
L
l
d
Ω
,
(20)
∫
ω
J
∩
ω
L
∇
Ψ
J
·
σ
J
k
T
·
∇
Ψ
L
·
g
L
k
d
Ω
,
(21)
∫
ω
J
∩
ω
L
∇
Ψ
J
·
σ
J
k
T
·
∇
Ψ
L
·
σ
L
l
d
Ω
,where g 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 [33].Theorem 4 (letΩ be a 2D convex polygon).
Also let{
Ψ
J
(
x
,
y
)
:
J
=
1,2
,
…
,
N
} be a partition of unity corresponding to the patch {
Q
J
:
J
=
1,2
,
…
,
N
} of Ω. Also denote Q
J
f
l
a
t for the flat-top region of Ψ
J
(
x
,
y
). Assume there exists a positive integer M such that, for every x
∈
Ω, card
⁡
{
J
:
x
∈
ω
J
}
≤
M, where ω
J
=
s
u
p
p
(
Ψ
J
). Let a collection of local approximation spaces V
J
⊂
H
1
(
Ω
∩
ω
J
) be given on each patch Q
J.
Let us assume the local approximation spaceV
J has the following approximation properties on each patch Q
J.
The functionu
∈
H
1
(
Ω
) can be approximated by a function g
J
∈
V
J such that (1)
u
-
g
J
L
2
(
Ω
∩
w
J
)
≤
ϵ
J
(
0
),
(2)
u
-
g
J
L
2
(
Ω
∩
(
ω
J
∖
ω
J
f
l
a
t
)
)
≤
ϵ
J
,
n
o
n
-
f
l
a
t
(
0
),
(3)
∇
(
u
-
g
J
)
L
2
(
Ω
∩
w
J
)
≤
ϵ
J
(
1
),for all J
=
1,2
,
…
,
N
. Then,(22)
u
∗
=
∑
J
=
1
N
Ψ
J
·
g
J satisfies the following global estimates:(23)
u
-
u
∗
L
2
Ω
2
≤
2
M
∑
J
=
1
N
ϵ
J
0
2
,
∇
u
-
u
∗
L
2
Ω
2
≤
4
M
∑
J
=
1
N
C
∇
δ
2
ϵ
J
,
n
o
n
-
f
l
a
t
0
2
+
ϵ
J
1
2
,where C
∇ is a constant independent of δ.


### 2.3. 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.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 width2
δ.Figure 3
Shaded rectangle is the enriched patch, which is subdivided by disjoint regionsD and D
c for numerical integration.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.Theorem 5.
Consider(24)
∫
Ω
flat
∇
x
,
y
g
i
r
,
θ
T
∇
x
,
y
g
j
r
,
θ
d
x
d
y
=
1
2
+
i
1
2
+
j
∫
Ω
flat
r
i
+
j
cos
⁡
i
-
j
θ
d
r
d
θ
,where Ω
flat denotes the flat-top region of the PU function corresponding to the enriched patch shown in Figure 3 and g
k
(
r
,
θ
)
=
r
k
+
1
/
2
cos
⁡
(
(
k
+
1
/
2
)
θ
).Proof.
By a coordinate transformation, the gradient operator is transformed as follows:(25)
∇
x
,
y
=
∂
∂
x
∂
∂
x
=
cos
⁡
θ
-
1
r
sin
⁡
θ
sin
⁡
θ
1
r
cos
⁡
θ
∇
r
,
θ
=
M
∂
∂
r
∂
∂
θ
,where M
=
cos
⁡
θ
-
(
1
/
r
)
sin
⁡
θ
sin
⁡
θ
(
1
/
r
)
cos
⁡
θ. For simplicity, we suppress (
r
,
θ
) in g
(
r
,
θ
). Then,(26)
∫
Ω
flat
∇
x
,
y
g
i
r
,
θ
T
∇
x
,
y
g
j
r
,
θ
d
x
d
y
=
∫
Ω
flat
M
∇
r
,
θ
g
i
T
M
∇
r
,
θ
g
j
r
d
r
d
θ
=
∫
Ω
flat
∇
r
,
θ
g
i
T
M
T
M
∇
r
,
θ
g
j
r
d
r
d
θ
=
∫
Ω
flat
∂
g
i
∂
r
,
∂
g
i
∂
θ
1
0
0
1
r
2
∂
g
j
∂
r
∂
g
j
∂
θ
r
d
r
d
θ
=
∫
Ω
flat
∂
g
i
∂
r
∂
g
j
∂
r
+
1
r
2
∂
g
i
∂
θ
∂
g
j
∂
θ
r
d
r
d
θ
.Since g
k
(
r
,
θ
)
=
r
1
/
2
+
k
cos
⁡
1
/
2
+
k
θ, (27)
∂
g
i
∂
r
=
1
2
+
i
r
i
-
1
/
2
cos
⁡
1
2
+
i
θ
,
∂
g
j
∂
r
=
1
2
+
j
r
j
-
1
/
2
cos
⁡
1
2
+
j
θ
,
∂
g
i
∂
θ
=
-
1
2
+
i
r
1
/
2
+
i
sin
⁡
1
2
+
i
θ
,
∂
g
j
∂
θ
=
-
1
2
+
j
r
1
/
2
+
j
sin
⁡
1
2
+
j
θ
.Plugging (27) into (26), we have the following:(28)
∫
Ω
flat
∇
x
,
y
g
i
r
,
θ
T
∇
x
,
y
g
j
r
,
θ
d
x
d
y
=
∫
Ω
flat
1
2
+
i
1
2
+
j
r
i
+
j
cos
⁡
1
2
+
i
θ
cos
⁡
1
2
+
j
θ
d
r
d
θ
+
∫
Ω
flat
1
2
+
i
1
2
+
j
r
i
+
j
sin
⁡
1
2
+
i
θ
sin
⁡
1
2
+
j
θ
d
r
d
θ
=
1
2
+
i
1
2
+
j
∫
Ω
flat
r
i
+
j
cos
⁡
1
2
+
i
θ
-
1
2
+
j
θ
d
r
d
θ
=
1
2
+
i
1
2
+
j
∫
Ω
flat
r
i
+
j
c
o
s
i
-
j
θ
d
r
d
θ
, which proves the theorem.LetI
,
J be the global indices corresponding to the local enrichment functions g
i
,
g
j, respectively, and let S
I
J be (
I
,
J
)
t
h component of stiffness matrix S. Then S
I
J is given by the following integration:(29)
S
I
J
=
∫
O
i
∩
O
j
∇
x
,
y
Ψ
i
x
,
y
g
i
r
,
θ
T
∇
x
,
y
Ψ
j
x
,
y
g
j
r
,
θ
d
x
d
y
,where O
i and O
j are the supports of the partition of unity functions Ψ
i and Ψ
j, respectively. Recall that the functions g
i and g
j are globally defined functions. Then, the overlapping region O
i
∩
O
j can be decomposed as O
flat
∪
˙
O
nonflat: a disjoint union of O
flat and O
nonflat. Since ∇
Ψ
i
=
∇
Ψ
j
=
(
0,0
)
T on O
flat, S
I
J of (29) can be written: (30)
S
I
J
=
∫
O
flat
∇
x
,
y
Ψ
i
g
i
T
∇
x
,
y
Ψ
j
g
j
d
x
d
y
+
∫
O
nonflat
∇
x
,
y
Ψ
i
g
i
T
∇
x
,
y
Ψ
j
g
j
d
x
d
y
=
∫
O
flat
∇
x
,
y
g
i
T
∇
x
,
y
g
j
d
x
d
y
︸
A
+
∫
O
nonflat
∇
x
,
y
Ψ
i
g
i
T
∇
x
,
y
Ψ
j
g
j
d
x
d
y
︸
B
,where Ψ
i
=
Ψ
i
(
x
,
y
)
,
Ψ
j
=
Ψ
j
(
x
,
y
) and g
i
=
g
i
(
r
,
θ
)
,
g
j
=
g
j
(
r
,
θ
).Let us examine(
B
) first. Note that Ψ
j is a piecewise polynomial and the integral region O
nonflat is far away from the singularity of g
j. Also it is important to note that the integral region has a rectangular shape. Therefore, integral (
B
) is straight forward to evaluate. With enough Gauss points, integral (
B
) should be nearly exact.LetP
:
(
r
,
θ
)
→
(
x
,
y
) be the coordinate transformation, x
=
r
cos
⁡
θ and
y
=
r
sin
⁡
θ, and let T
R
:
(
s
,
t
)
→
(
x
,
y
) be the patch mapping from the reference patch P to a quadrangle patch R. See Figure 4(a). Then we have(31)
∂
∂
x
Ψ
i
g
i
=
∂
Ψ
i
∂
x
g
i
+
Ψ
i
∂
g
i
∂
r
∂
r
∂
x
+
∂
g
i
∂
θ
∂
θ
∂
x
=
∂
Ψ
i
∂
x
g
i
∘
P
-
1
+
Ψ
i
∂
g
i
∂
r
cos
⁡
θ
-
∂
g
i
∂
θ
1
r
sin
⁡
θ
∘
P
-
1
.Similarly,(32)
∂
∂
x
Ψ
j
g
j
=
∂
Ψ
j
∂
x
g
j
∘
P
-
1
+
Ψ
j
∂
g
j
∂
r
cos
⁡
θ
-
∂
g
j
∂
θ
1
r
sin
⁡
θ
∘
P
-
1
,
∂
∂
y
Ψ
i
g
i
=
∂
Ψ
i
∂
y
g
i
∘
P
-
1
+
Ψ
i
∂
g
i
∂
r
sin
⁡
θ
+
∂
g
i
∂
θ
1
r
cos
⁡
θ
∘
P
-
1
,
∂
∂
y
Ψ
j
g
j
=
∂
Ψ
j
∂
y
g
j
∘
P
-
1
+
Ψ
j
∂
g
j
∂
r
sin
⁡
θ
+
∂
g
j
∂
θ
1
r
cos
⁡
θ
∘
P
-
1
. Note ∂
/
∂
x
Ψ
j
g
j, ∂
/
∂
x
Ψ
j
g
j, ∂
/
∂
y
Ψ
i
g
i, and ∂
/
∂
y
Ψ
j
g
j are all function of x
,
y. Let us define G
(
x
,
y
)
=
∂
/
∂
x
Ψ
i
g
i
∂
/
∂
x
Ψ
j
g
j
+
∂
/
∂
y
Ψ
i
g
i
∂
/
∂
y
Ψ
j
g
j. Then, (33)
B
=
∫
O
nonflat
∂
∂
x
Ψ
i
g
i
∂
∂
x
Ψ
j
g
j
+
∂
∂
y
Ψ
i
g
i
∂
∂
y
Ψ
j
g
j
d
x
d
y
=
∫
O
nonflat
G
x
,
y
d
x
d
y
=
∫
P
G
∘
T
R
J
T
R
d
s
d
t
=
∑
n
=
1
N
g
∑
m
=
1
N
g
G
∘
T
B
k
g
n
,
g
m
J
T
B
R
g
n
,
g
m
W
n
W
m
,where N
g is the number of Legendre-Gauss quadrature points, g
n is the Legendre-Gauss point, and W
n is the nth quadrature weight.Figure 4
Two patch mappings.
(a)
A rectangular patch mappingT
R
(b)
A blending mappingT
BOn the other hand, by Theorem5, Part (
A
) is as follows: (34)
A
=
1
2
+
i
1
2
+
j
∫
O
flat
r
i
+
j
cos
⁡
i
-
j
θ
d
r
d
θ
.
O
flat of the enriched patch can be further decomposed into two separate regions D and D
c. D
=
{
(
r
,
θ
)
∣
0
≤
r
≤
R
,
0
≤
θ
≤
π
} in Figure 3 is the half disk and D
c is its complement; that is, O
flat
=
D
∪
˙
D
c (see Figure 3). On D, the integral of part (
A
) can be obtained exactly.Then, integral(
A
) over D can be simplified as follows:(35)
1
2
+
i
π
R
2
i
+
1
2
if
i
=
j
,
0
if
i
≠
j
.Using (35), we have (36)
A
=
1
2
+
i
π
R
2
i
+
1
2
δ
j
i
+
1
2
+
i
1
2
+
j
∫
D
c
r
i
+
j
cos
⁡
i
-
j
θ
d
r
d
θ
︸
C
,where δ
j
i is the Kronecker delta. Now over D
c region, we need to perform numerical integration because of its irregular shape. Since D
c region is arch-shaped, we need to use the patch mapping of the blending type for this integration.Remark 6.
Note that the integrandcos
⁡
i
-
j
θ becomes highly oscillatory as the value of i
-
j increases. For instance, for i
=
40
,
j
=
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 mappingT
B
:
(
s
,
t
)
→
(
x
,
y
), Figure 4(b), maps the reference patch P onto a quadrangle Q with a curved side. Let P
1
=
(
R
,
θ
2
) and P
2
=
(
R
,
θ
1
) in terms of polar coordinates and P
3
=
(
x
3
,
y
3
) and P
4
=
(
x
4
,
y
4
) in rectangular coordinates. The blending patch mapping T
B is defined as follows: (37)
x
=
1
-
t
2
R
cos
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
+
1
+
t
2
x
3
1
+
s
2
+
x
4
1
-
s
2
,
y
=
1
-
t
2
R
sin
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
+
1
+
t
2
y
3
1
+
s
2
+
y
4
1
-
s
2
.Also its partial derivatives are as follows: (38)
∂
x
∂
s
=
t
-
1
2
R
sin
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
θ
1
-
θ
2
2
+
1
+
t
4
x
3
-
x
4
,
(39)
∂
x
∂
t
=
-
1
2
R
cos
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
+
1
2
x
3
1
+
s
2
+
x
4
1
-
s
2
,
(40)
∂
y
∂
s
=
1
-
t
2
R
cos
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
θ
1
-
θ
2
2
+
1
+
t
4
y
3
-
y
4
,
(41)
∂
y
∂
t
=
-
1
2
R
sin
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
+
1
2
y
3
1
+
s
2
+
y
4
1
-
s
2
. The determinant of Jacobian (T
B), denoted by J
T
B, is ∂
x
/
∂
s
∂
y
/
∂
t
-
∂
x
/
∂
t
∂
y
/
∂
s. Hence integral over region Q can be transformed over region P: (42)
∫
Q
F
x
,
y
d
x
d
y
=
∫
P
F
∘
T
B
∂
x
∂
s
∂
y
∂
t
-
∂
x
∂
t
∂
y
∂
s
d
s
d
t
. As pointed out earlier to obtain accurate numerical integration, D
c region should be partitioned into many pieces to capture the oscillatory behavior of the integrand. Let us say D
c is partitioned into N pieces. Also let T
B
k be the transformation from the reference square P into the integral region D
k
c, where D
c
=
⋃
k
=
1
N
D
k
c. Then integral (
C
) becomes the following: (43)
C
=
∑
k
=
1
N
∫
D
k
c
r
i
+
j
cos
⁡
i
-
j
θ
d
r
d
θ
=
∑
k
=
1
N
∫
P
r
i
+
j
cos
⁡
i
-
j
θ
∘
T
B
k
J
T
B
k
d
s
d
t
=
∑
k
=
1
N
∬
-
1
1
r
i
+
j
cos
⁡
i
-
j
θ
∘
T
B
k
J
T
B
k
d
s
d
t
=
∑
k
=
1
N
∑
n
=
1
N
g
∑
m
=
1
N
g
r
i
+
j
cos
⁡
i
-
j
θ
∘
T
B
k
g
n
,
g
m
J
T
B
k
g
n
,
g
m
W
n
W
m
,where N
g denotes the number of Legendre-Gauss quadrature points, g
n denotes the Legendre-Gauss point, and W
n is the nth quadrature weight. Thus, we have completed the calculation of S
I
J. 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.


## 2.1. 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.Definition 2.
The enriched meshfree approximation spaceV
enrich is defined as follows. Suppose one wants to enrich local approximation functions on the patches (11)
Q
J
1
,
Q
J
2
,
…
,
Q
J
swith sets of special functions: (12)
σ
J
i
,
k
,
k
=
1,2
,
…
,
N
J
i
for
i
=
1,2
,
…
,
s
.Then V
enrich is the vector space spanned by (13)
V
mfree
∪
Ψ
J
i
σ
J
i
k
:
k
=
1,2
,
…
,
N
J
i
,
i
=
1,2
,
…
,
s
.The enriched meshfree method for the second-order elliptic equation with the boundary condition (1) associated with the enriched meshfree approximation space V
enrich is as follows.Find(14)
u
enrich
=
∑
J
=
1
M
∑
k
=
1
N
J
c
J
k
Ψ
J
g
J
k
+
∑
i
=
1
s
∑
l
=
1
N
J
i
d
J
i
,
l
Ψ
J
i
σ
J
i
,
l
,such that(15)
B
u
enrich
,
v
=
F
v
,
for
v
∈
V
enrich
,where B and F are defined as (4).Calculation of a stiffness matrix involves the following integration:(16)
∫
ω
J
∩
ω
L
∇
v
J
,
i
T
·
∇
v
L
,
j
d
Ω
,where J and L indicate the patch number and ∇
=
(
∂
/
∂
x
,
∂
/
∂
y
)
T. Also v
J
,
i
,
v
L
,
j
∈
V
enrich and the components of the load vector are given as follows:(17)
∫
ω
J
∩
ω
L
f
·
v
J
,
k
d
Ω
.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:(18)
∫
ω
J
∩
ω
L
∇
Ψ
J
·
g
J
k
T
·
∇
Ψ
L
·
g
L
l
d
Ω
,
(19)
∫
ω
J
∩
ω
L
∇
Ψ
J
·
g
J
k
T
·
∇
Ψ
L
·
σ
L
l
d
Ω
,
(20)
∫
ω
J
∩
ω
L
∇
Ψ
J
·
σ
J
k
T
·
∇
Ψ
L
·
g
L
k
d
Ω
,
(21)
∫
ω
J
∩
ω
L
∇
Ψ
J
·
σ
J
k
T
·
∇
Ψ
L
·
σ
L
l
d
Ω
,where g 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 [33].Theorem 4 (letΩ be a 2D convex polygon).
Also let{
Ψ
J
(
x
,
y
)
:
J
=
1,2
,
…
,
N
} be a partition of unity corresponding to the patch {
Q
J
:
J
=
1,2
,
…
,
N
} of Ω. Also denote Q
J
f
l
a
t for the flat-top region of Ψ
J
(
x
,
y
). Assume there exists a positive integer M such that, for every x
∈
Ω, card
⁡
{
J
:
x
∈
ω
J
}
≤
M, where ω
J
=
s
u
p
p
(
Ψ
J
). Let a collection of local approximation spaces V
J
⊂
H
1
(
Ω
∩
ω
J
) be given on each patch Q
J.
Let us assume the local approximation spaceV
J has the following approximation properties on each patch Q
J.
The functionu
∈
H
1
(
Ω
) can be approximated by a function g
J
∈
V
J such that (1)
u
-
g
J
L
2
(
Ω
∩
w
J
)
≤
ϵ
J
(
0
),
(2)
u
-
g
J
L
2
(
Ω
∩
(
ω
J
∖
ω
J
f
l
a
t
)
)
≤
ϵ
J
,
n
o
n
-
f
l
a
t
(
0
),
(3)
∇
(
u
-
g
J
)
L
2
(
Ω
∩
w
J
)
≤
ϵ
J
(
1
),for all J
=
1,2
,
…
,
N
. Then,(22)
u
∗
=
∑
J
=
1
N
Ψ
J
·
g
J satisfies the following global estimates:(23)
u
-
u
∗
L
2
Ω
2
≤
2
M
∑
J
=
1
N
ϵ
J
0
2
,
∇
u
-
u
∗
L
2
Ω
2
≤
4
M
∑
J
=
1
N
C
∇
δ
2
ϵ
J
,
n
o
n
-
f
l
a
t
0
2
+
ϵ
J
1
2
,where C
∇ is a constant independent of δ.


## 2.3. 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.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 width2
δ.Figure 3
Shaded rectangle is the enriched patch, which is subdivided by disjoint regionsD and D
c for numerical integration.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.Theorem 5.
Consider(24)
∫
Ω
flat
∇
x
,
y
g
i
r
,
θ
T
∇
x
,
y
g
j
r
,
θ
d
x
d
y
=
1
2
+
i
1
2
+
j
∫
Ω
flat
r
i
+
j
cos
⁡
i
-
j
θ
d
r
d
θ
,where Ω
flat denotes the flat-top region of the PU function corresponding to the enriched patch shown in Figure 3 and g
k
(
r
,
θ
)
=
r
k
+
1
/
2
cos
⁡
(
(
k
+
1
/
2
)
θ
).Proof.
By a coordinate transformation, the gradient operator is transformed as follows:(25)
∇
x
,
y
=
∂
∂
x
∂
∂
x
=
cos
⁡
θ
-
1
r
sin
⁡
θ
sin
⁡
θ
1
r
cos
⁡
θ
∇
r
,
θ
=
M
∂
∂
r
∂
∂
θ
,where M
=
cos
⁡
θ
-
(
1
/
r
)
sin
⁡
θ
sin
⁡
θ
(
1
/
r
)
cos
⁡
θ. For simplicity, we suppress (
r
,
θ
) in g
(
r
,
θ
). Then,(26)
∫
Ω
flat
∇
x
,
y
g
i
r
,
θ
T
∇
x
,
y
g
j
r
,
θ
d
x
d
y
=
∫
Ω
flat
M
∇
r
,
θ
g
i
T
M
∇
r
,
θ
g
j
r
d
r
d
θ
=
∫
Ω
flat
∇
r
,
θ
g
i
T
M
T
M
∇
r
,
θ
g
j
r
d
r
d
θ
=
∫
Ω
flat
∂
g
i
∂
r
,
∂
g
i
∂
θ
1
0
0
1
r
2
∂
g
j
∂
r
∂
g
j
∂
θ
r
d
r
d
θ
=
∫
Ω
flat
∂
g
i
∂
r
∂
g
j
∂
r
+
1
r
2
∂
g
i
∂
θ
∂
g
j
∂
θ
r
d
r
d
θ
.Since g
k
(
r
,
θ
)
=
r
1
/
2
+
k
cos
⁡
1
/
2
+
k
θ, (27)
∂
g
i
∂
r
=
1
2
+
i
r
i
-
1
/
2
cos
⁡
1
2
+
i
θ
,
∂
g
j
∂
r
=
1
2
+
j
r
j
-
1
/
2
cos
⁡
1
2
+
j
θ
,
∂
g
i
∂
θ
=
-
1
2
+
i
r
1
/
2
+
i
sin
⁡
1
2
+
i
θ
,
∂
g
j
∂
θ
=
-
1
2
+
j
r
1
/
2
+
j
sin
⁡
1
2
+
j
θ
.Plugging (27) into (26), we have the following:(28)
∫
Ω
flat
∇
x
,
y
g
i
r
,
θ
T
∇
x
,
y
g
j
r
,
θ
d
x
d
y
=
∫
Ω
flat
1
2
+
i
1
2
+
j
r
i
+
j
cos
⁡
1
2
+
i
θ
cos
⁡
1
2
+
j
θ
d
r
d
θ
+
∫
Ω
flat
1
2
+
i
1
2
+
j
r
i
+
j
sin
⁡
1
2
+
i
θ
sin
⁡
1
2
+
j
θ
d
r
d
θ
=
1
2
+
i
1
2
+
j
∫
Ω
flat
r
i
+
j
cos
⁡
1
2
+
i
θ
-
1
2
+
j
θ
d
r
d
θ
=
1
2
+
i
1
2
+
j
∫
Ω
flat
r
i
+
j
c
o
s
i
-
j
θ
d
r
d
θ
, which proves the theorem.LetI
,
J be the global indices corresponding to the local enrichment functions g
i
,
g
j, respectively, and let S
I
J be (
I
,
J
)
t
h component of stiffness matrix S. Then S
I
J is given by the following integration:(29)
S
I
J
=
∫
O
i
∩
O
j
∇
x
,
y
Ψ
i
x
,
y
g
i
r
,
θ
T
∇
x
,
y
Ψ
j
x
,
y
g
j
r
,
θ
d
x
d
y
,where O
i and O
j are the supports of the partition of unity functions Ψ
i and Ψ
j, respectively. Recall that the functions g
i and g
j are globally defined functions. Then, the overlapping region O
i
∩
O
j can be decomposed as O
flat
∪
˙
O
nonflat: a disjoint union of O
flat and O
nonflat. Since ∇
Ψ
i
=
∇
Ψ
j
=
(
0,0
)
T on O
flat, S
I
J of (29) can be written: (30)
S
I
J
=
∫
O
flat
∇
x
,
y
Ψ
i
g
i
T
∇
x
,
y
Ψ
j
g
j
d
x
d
y
+
∫
O
nonflat
∇
x
,
y
Ψ
i
g
i
T
∇
x
,
y
Ψ
j
g
j
d
x
d
y
=
∫
O
flat
∇
x
,
y
g
i
T
∇
x
,
y
g
j
d
x
d
y
︸
A
+
∫
O
nonflat
∇
x
,
y
Ψ
i
g
i
T
∇
x
,
y
Ψ
j
g
j
d
x
d
y
︸
B
,where Ψ
i
=
Ψ
i
(
x
,
y
)
,
Ψ
j
=
Ψ
j
(
x
,
y
) and g
i
=
g
i
(
r
,
θ
)
,
g
j
=
g
j
(
r
,
θ
).Let us examine(
B
) first. Note that Ψ
j is a piecewise polynomial and the integral region O
nonflat is far away from the singularity of g
j. Also it is important to note that the integral region has a rectangular shape. Therefore, integral (
B
) is straight forward to evaluate. With enough Gauss points, integral (
B
) should be nearly exact.LetP
:
(
r
,
θ
)
→
(
x
,
y
) be the coordinate transformation, x
=
r
cos
⁡
θ and
y
=
r
sin
⁡
θ, and let T
R
:
(
s
,
t
)
→
(
x
,
y
) be the patch mapping from the reference patch P to a quadrangle patch R. See Figure 4(a). Then we have(31)
∂
∂
x
Ψ
i
g
i
=
∂
Ψ
i
∂
x
g
i
+
Ψ
i
∂
g
i
∂
r
∂
r
∂
x
+
∂
g
i
∂
θ
∂
θ
∂
x
=
∂
Ψ
i
∂
x
g
i
∘
P
-
1
+
Ψ
i
∂
g
i
∂
r
cos
⁡
θ
-
∂
g
i
∂
θ
1
r
sin
⁡
θ
∘
P
-
1
.Similarly,(32)
∂
∂
x
Ψ
j
g
j
=
∂
Ψ
j
∂
x
g
j
∘
P
-
1
+
Ψ
j
∂
g
j
∂
r
cos
⁡
θ
-
∂
g
j
∂
θ
1
r
sin
⁡
θ
∘
P
-
1
,
∂
∂
y
Ψ
i
g
i
=
∂
Ψ
i
∂
y
g
i
∘
P
-
1
+
Ψ
i
∂
g
i
∂
r
sin
⁡
θ
+
∂
g
i
∂
θ
1
r
cos
⁡
θ
∘
P
-
1
,
∂
∂
y
Ψ
j
g
j
=
∂
Ψ
j
∂
y
g
j
∘
P
-
1
+
Ψ
j
∂
g
j
∂
r
sin
⁡
θ
+
∂
g
j
∂
θ
1
r
cos
⁡
θ
∘
P
-
1
. Note ∂
/
∂
x
Ψ
j
g
j, ∂
/
∂
x
Ψ
j
g
j, ∂
/
∂
y
Ψ
i
g
i, and ∂
/
∂
y
Ψ
j
g
j are all function of x
,
y. Let us define G
(
x
,
y
)
=
∂
/
∂
x
Ψ
i
g
i
∂
/
∂
x
Ψ
j
g
j
+
∂
/
∂
y
Ψ
i
g
i
∂
/
∂
y
Ψ
j
g
j. Then, (33)
B
=
∫
O
nonflat
∂
∂
x
Ψ
i
g
i
∂
∂
x
Ψ
j
g
j
+
∂
∂
y
Ψ
i
g
i
∂
∂
y
Ψ
j
g
j
d
x
d
y
=
∫
O
nonflat
G
x
,
y
d
x
d
y
=
∫
P
G
∘
T
R
J
T
R
d
s
d
t
=
∑
n
=
1
N
g
∑
m
=
1
N
g
G
∘
T
B
k
g
n
,
g
m
J
T
B
R
g
n
,
g
m
W
n
W
m
,where N
g is the number of Legendre-Gauss quadrature points, g
n is the Legendre-Gauss point, and W
n is the nth quadrature weight.Figure 4
Two patch mappings.
(a)
A rectangular patch mappingT
R
(b)
A blending mappingT
BOn the other hand, by Theorem5, Part (
A
) is as follows: (34)
A
=
1
2
+
i
1
2
+
j
∫
O
flat
r
i
+
j
cos
⁡
i
-
j
θ
d
r
d
θ
.
O
flat of the enriched patch can be further decomposed into two separate regions D and D
c. D
=
{
(
r
,
θ
)
∣
0
≤
r
≤
R
,
0
≤
θ
≤
π
} in Figure 3 is the half disk and D
c is its complement; that is, O
flat
=
D
∪
˙
D
c (see Figure 3). On D, the integral of part (
A
) can be obtained exactly.Then, integral(
A
) over D can be simplified as follows:(35)
1
2
+
i
π
R
2
i
+
1
2
if
i
=
j
,
0
if
i
≠
j
.Using (35), we have (36)
A
=
1
2
+
i
π
R
2
i
+
1
2
δ
j
i
+
1
2
+
i
1
2
+
j
∫
D
c
r
i
+
j
cos
⁡
i
-
j
θ
d
r
d
θ
︸
C
,where δ
j
i is the Kronecker delta. Now over D
c region, we need to perform numerical integration because of its irregular shape. Since D
c region is arch-shaped, we need to use the patch mapping of the blending type for this integration.Remark 6.
Note that the integrandcos
⁡
i
-
j
θ becomes highly oscillatory as the value of i
-
j increases. For instance, for i
=
40
,
j
=
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 mappingT
B
:
(
s
,
t
)
→
(
x
,
y
), Figure 4(b), maps the reference patch P onto a quadrangle Q with a curved side. Let P
1
=
(
R
,
θ
2
) and P
2
=
(
R
,
θ
1
) in terms of polar coordinates and P
3
=
(
x
3
,
y
3
) and P
4
=
(
x
4
,
y
4
) in rectangular coordinates. The blending patch mapping T
B is defined as follows: (37)
x
=
1
-
t
2
R
cos
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
+
1
+
t
2
x
3
1
+
s
2
+
x
4
1
-
s
2
,
y
=
1
-
t
2
R
sin
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
+
1
+
t
2
y
3
1
+
s
2
+
y
4
1
-
s
2
.Also its partial derivatives are as follows: (38)
∂
x
∂
s
=
t
-
1
2
R
sin
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
θ
1
-
θ
2
2
+
1
+
t
4
x
3
-
x
4
,
(39)
∂
x
∂
t
=
-
1
2
R
cos
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
+
1
2
x
3
1
+
s
2
+
x
4
1
-
s
2
,
(40)
∂
y
∂
s
=
1
-
t
2
R
cos
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
θ
1
-
θ
2
2
+
1
+
t
4
y
3
-
y
4
,
(41)
∂
y
∂
t
=
-
1
2
R
sin
⁡
θ
1
-
θ
2
2
s
+
1
+
θ
2
+
1
2
y
3
1
+
s
2
+
y
4
1
-
s
2
. The determinant of Jacobian (T
B), denoted by J
T
B, is ∂
x
/
∂
s
∂
y
/
∂
t
-
∂
x
/
∂
t
∂
y
/
∂
s. Hence integral over region Q can be transformed over region P: (42)
∫
Q
F
x
,
y
d
x
d
y
=
∫
P
F
∘
T
B
∂
x
∂
s
∂
y
∂
t
-
∂
x
∂
t
∂
y
∂
s
d
s
d
t
. As pointed out earlier to obtain accurate numerical integration, D
c region should be partitioned into many pieces to capture the oscillatory behavior of the integrand. Let us say D
c is partitioned into N pieces. Also let T
B
k be the transformation from the reference square P into the integral region D
k
c, where D
c
=
⋃
k
=
1
N
D
k
c. Then integral (
C
) becomes the following: (43)
C
=
∑
k
=
1
N
∫
D
k
c
r
i
+
j
cos
⁡
i
-
j
θ
d
r
d
θ
=
∑
k
=
1
N
∫
P
r
i
+
j
cos
⁡
i
-
j
θ
∘
T
B
k
J
T
B
k
d
s
d
t
=
∑
k
=
1
N
∬
-
1
1
r
i
+
j
cos
⁡
i
-
j
θ
∘
T
B
k
J
T
B
k
d
s
d
t
=
∑
k
=
1
N
∑
n
=
1
N
g
∑
m
=
1
N
g
r
i
+
j
cos
⁡
i
-
j
θ
∘
T
B
k
g
n
,
g
m
J
T
B
k
g
n
,
g
m
W
n
W
m
,where N
g denotes the number of Legendre-Gauss quadrature points, g
n denotes the Legendre-Gauss point, and W
n is the nth quadrature weight. Thus, we have completed the calculation of S
I
J. 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.Table1 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 u
∈
H
1
(
Ω
) is U
=
1
/
2
B
(
u
,
u
). The relative error in the energy norm is defined by (44)
U
∗
-
U
Li
U
Li
1
/
2
,where U
∗ and U
L
i are the strain energy obtained by the solution of proposed method and provided in the literature [9], respectively.Table 1
Enriched meshfree method errors, measured in max and energy norm.
p degree
dof
err
⁡
m
ax
err
⁡
energy
2
349
2.75
×
1
0
-
3
7.24
×
1
0
-
5
4
940
5.69
×
1
0
-
7
1.87
×
1
0
-
6
6
1830
2.13
×
1
0
-
8
1.74
×
1
0
-
7
8
3019
5.37
×
1
0
-
9
4.37
×
1
0
-
8Figure5 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.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.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–36] in the framework of the p-version FEM. Table 2 compares the result of the proposed method with RSPM and p-FEM with MAM. We conclude that the enriched meshfree method is superior to both RSPM and p-FEM with MAM for the Motz problem. Also it is important to note that both p-version FEM with MAM and RSPM can not estimate d
~
i in (3) but the proposed method can.Table 2
Relative errors are compared for differentp degrees in energy norm.
p degree
p-FEM withMAM
MeshfreeRSPM
Proposedmethod
2
9.00
×
1
0
-
3
1.54
×
1
0
-
2
7.24
×
1
0
-
5
4
2.38
×
1
0
-
4
3.17
×
1
0
-
4
1.87
×
1
0
-
6
6
3.42
×
1
0
-
6
2.40
×
1
0
-
5
1.74
×
1
0
-
7
8
4.83
×
1
0
-
8
7.23
×
1
0
-
7
4.37
×
1
0
-
8We 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 Table4, 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 d
i in (14) are not intended to be used in (3). However, d
~
i’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 Ψ
J
i, in (14), has the value of one if we restrict the numerical solution locally near the singularity. Therefore, we may compare d
i in (14) with Li’s solution [9], which is the best known numerical solution for the Motz problem. We denote the coefficients from our method to d
~
i
∗ and Li’s solution to d
~
i
†. d
~
i, i
=
1
,
…
,
40, is listed in Table 3. This time we compare the errors of ∑
i
=
0
39
d
~
i
∗
r
i
+
1
/
2
cos
⁡
⁡
(
i
+
1
/
2
)
θ and ∑
i
=
0
34
d
~
i
†
r
i
+
1
/
2
cos
⁡
(
i
+
1
/
2
)
θ along the boundary Γ
2 where the true solution is 500. As mentioned earlier, the enriched meshfree solution given in (14) does not have any error on the boundary Γ
2. To compare the accuracy with the best known solution, d
~
i
∗,  i
=
1
,
…
,
40, are extracted from the meshfree solution and applied to (3). Therefore, the series solution obtained by using d
~
i, i
=
1
,
…
,
40, as coefficients exhibits error on Γ
2 and they are measured by 500
-
∑
i
=
0
39
d
~
i
∗
r
i
+
1
/
2
cos
⁡
i
+
1
/
2
θ and 500
-
∑
i
=
0
34
d
~
i
†
r
i
+
1
/
2
cos
⁡
i
+
1
/
2
θ, respectively. The errors along the boundary Γ
2 are plotted in Figure 6.Table 3
Computed coefficientsd
~
i
∗ of ∑
i
=
1
40
d
~
i
∗
r
i
-
1
/
2
c
o
s
(
(
i
-
1
/
2
)
θ
) with the proposed enriched meshfree method.
ith
Sign
Coefficients
1
+
401.1624537452345
2
+
87.65592019508790
3
+
17.23791507944674
4
-
8.071215259682478
5
+
1.440272717022879
6
+
3.310548859207677
e
-
01
7
+
2.754373445087275
e
-
01
8
-
8.693299465893235
e
-
02
9
+
3.360487861931393
e
-
02
10
+
1.538437428404088
e
-
02
11
+
7.302301728215649
e
-
03
12
-
3.184113774923665
e
-
03
13
+
1.220646385725587
e
-
03
14
+
5.309652968145126
e
-
04
15
+
2.715121558566150
e
-
04
16
-
1.200455757406748
e
-
04
17
+
5.053973008462738
e
-
05
18
+
2.316701639900009
e
-
05
19
+
1.153538904754159
e
-
05
20
-
5.293975100353315
e
-
06
21
+
2.290744823584897
e
-
06
22
+
1.063692985558411
e
-
06
23
+
5.314980562968453
e
-
07
24
-
2.454035587720472
e
-
07
25
+
1.094258200815170
e
-
07
26
+
5.204123234110838
e
-
08
27
+
2.594548101382485
e
-
08
28
-
1.091753910227160
e
-
08
29
+
5.274360833540486
e
-
09
30
+
2.733609395521398
e
-
09
31
+
1.354076917219339
e
-
09
32
-
1.816851682230465
e
-
10
33
+
2.266480380995557
e
-
10
34
+
1.541467769079195
e
-
10
35
+
7.396367522015557
e
-
11
36
+
5.058240306963314
e
-
11
37
+
6.165507191581204
e
-
12
38
+
6.544346229722119
e
-
12
39
+
2.943942531164191
e
-
12
40
+
5.127853852958069
e
-
12Table 4
Comparison of the first 25 coefficients obtained by other methods.
Exact
Li et al.
ISBFM
SFBIM
Trefftz
Enriched meshfree
i
N
=
20
N
=
35
N
=
40
N
=
75
N
=
35
N
=
40
[2]
[9]
[31]
[4]
[12]
Present
1
401.1624537452
401.1624537450
401.162453745
401.1624537452
401.1624537452344
401.1624537452345
2
87.6559201951
87.6559201941
87.6559202595
87.6559201951
87.65592019508793
87.65592019508790
3
17.2379150794
17.2379150819
17.2379150363
17.2379150794
17.23791507944678
17.23791507944674
4
-
8.0712152597
-
8.0712152607
-
8.0712152597
-
8.0712152597
-
8.071215259698
-
8.071215259682
5
1.4402727170
1.4402727163
1.4402727171
1.4402727170
1.440272717022
1.440272717022
6
0.3310548859
0.3310548866
0.3310548859
0.3310548859
0.33105488592
0.33105488592
7
0.2754373445
0.2754373447
0.2754373443
0.2754373445
0.27543734450
0.27543734450
8
-
0.0869329945
-
0.0869329948
-
0.0869329946
-
0.0869329945
-
0.08693299450
-
0.08693299465
9
0.0336048784
0.0336048781
-
0.0336048784
0.0336048784
0.03360487840
0.03360487861
10
0.0153843745
0.0153843747
0.0153843745
0.0153843745
0.01538437445
0.15384374284
11
0.0073023017
0.0073023019
0.0073023017
0.0073023017
0.00730230164
0.00730230172
12
-
0.0031841139
-
0.0031841138
-
0.0031841139
-
0.0031841139
-
0.00318411361
-
0.00318411377
13
0.0012206461
0.0012206456
0.0012206461
0.0012206461
0.00122064586
0.00122064638
14
0.0005309655
0.0005309655
0.0005309655
0.0005309655
0.00053096529
0.00053096529
15
0.0002715122
0.0002715122
0.0002715122
0.0002715122
0.00027151202
0.00027151215
16
-
0.0001200463
-
0.0001200450
-
0.0001200464
-
0.0001200464
-
0.00012004504
-
0.00012004557
17
0.0000505400
0.0000505387
0.0000505398
0.0000505398
0.00005053892
0.00005053973
18
0.000023167
0.0000231664
0.0000231668
0.0000231669
0.00002316625
0.00002316701
19
0.000011535
0.0000115349
0.0000115352
0.0000115353
0.00001153484
0.00001153538
20
-
0.000005295
-
0.0000052931
-
0.0000052957
-
0.0000052958
-
0.00000529323
-
0.00000529397
21
0.0000022895
0.0000022911
0.0000022911
0.00000228975
0.00000229744
22
0.0000010624
0.0000010632
0.0000010635
0.00000106239
0.00000106369
23
0.0000005307
0.0000005312
0.0000005314
0.00000053072
0.00000053149
24
-
0.0000002449
-
0.0000002473
-
0.0000002474
-
0.00000024507
-
0.00000024540
25
0.0000001085
0.0000001097
0.0000001087
0.00000010864
0.00000010942Figure 6
Errors alongΓ
2, where the Dirichlet boundary condition u
=
500 is specified.Remark 8.
The absolute error in maximum norm onΓ
2 is measured as 6.89
×
1
0
-
9 for the series solution with d
~
i
∗ and 5.47
×
1
0
-
9 for the series solution with d
~
i
†. In the mean square norm, the errors are 2.78
×
1
0
-
9 and 2.26
×
1
0
-
9, respectively. Although Li’s coefficients result in series solution that is 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 coefficientsd
~
j 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
~
i. 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 typer
α with α
<
1) but also oscillating singularities (of type r
α
sin
⁡
(
ϵ
log
⁡
r
) [36]).



---

*Source: 6324754-2016-07-10.xml*