# A New Iterative Method for Linear Systems from XFEM

Authors: Jianping Zhao; Yanren Hou; Haibiao Zheng; Bo Tang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The extend finite element method (XFEM) is popular in structural mechanics when dealing with the problem of the cracked domains. XFEM ends up with a linear system. However, XFEM usually leads to nonsymmetric and ill-conditioned stiff matrix. In this paper, we take the linear elastostatics governing equations as the model problem.
We propose a new iterative method to solve the linear equations. Here we separate two variablesU and Enr, so that we change the problems into solving the smaller scale equations iteratively. The new program can be easily applied. Finally, numerical examples show
that the proposed method is more efficient than common methods; we compare the L
2-error and the CPU time in whole process. Furthermore, the new XFEM can be applied and optimized in many other problems.


---

## Body

## 1. Introduction

Many physical phenomena can be modeled by partial differential equations with singularities and interfaces. Much attention has been focused on the development of interface problem in fluid dynamics and material science. The standard finite difference and finite element methods may not be successful in giving satisfactory numerical results for such problems. Hence, many new methods have been developed. Some of them are developed with the modifications in the standard methods, so that they can deal with the discontinuities and the singularities. Peskin developed immersed boundary method (IBM) [1] basically to model blood flow in the heart. Immersed interface method (IIM) [2–6] was basically developed because of the need for a second-order accurate finite difference method to solve linear elliptic, parabolic partial differential equations in which an integral equation may not be available. IIM will give a nonsymmetric system of equations. The extended finite element method (XFEM) [7, 8] enables the accurate approximation of solutions that involve jumps, kinks, singularities, and other locally nonsmooth features within elements. This is achieved by enriching the polynomial approximation space of the classical finite element method. The generalized finite element method (GEFM)/XFEM [9–11] has shown its potential in a variety of applications that involve nonsmooth solutions near interfaces. Among them are the simulation of cracks, shear bands, dislocations, solidification, and multifield problems. A large number of examples can be found in the real world where field quantities change rapidly over length scales that are small with respect to the observed domain. For the modeling of such phenomena, the resulting solutions typically involve discontinuities, singularities, high gradients, or other nonsmooth properties. For the approximation of nonsmooth solutions, one of the approaches is to enrich a polynomial approximation space such that the nonsmooth solutions can be modeled independently of the mesh. Methods that extend polynomial approximation spaces are called “enriched methods” in this work. The enrichment can be achieved by adding special shape functions (which are customized so as to capture jumps, singularities, etc.) to the polynomial approximation space. The total stiffness matrix consists of four parts: the FEM part, the enrich nodes part, and another two parts are the connecting parts between the real nodes and enrich nodes. Compared with FEM, the matrix size increases O
(
N
), while the FEM stiff matrix size is O
(
N
2
)
(
N
~
O
(
h
-
1
)
). But engineering calculations in XFEM always encountered the problem of solving ill-conditioned matrix. It made the XFEM widely.In this paper, according to the characteristics of the stiffness matrix, we use an iterative method to avoid dealing with complicated stiff matrix. First we can use the result of linear equation from standard FEM. The result can be used to be the initial value for iteration. Then we decompose the unknown values to two parts; each part is governed by two linear equations, whose dimension is smaller than the whole stiff matrix. The new programm can be easily applied. Finally, numerical examples show that the proposed method is more efficient than common methods. Furthermore, the new method works well for some problems while the common method (such as LU) fails. The new iterative method can be applied to many other problems.The outline of this paper is as follows. Section2 introduces the XFEM. Subsequently, the main part of this paper is Section 3 that we describe a new efficient and feasible method. We also analyze the time complexity of the algorithm. In Section 4, we briefly mention details on numerical examples. A final conclusion is drawn in Section 5.


## 2. The Extended Finite Element Method

The basic mathematical foundation of the partition of unity finite element method (PUFEM) was discussed in [12, 13]. They illustrated that PUFEM can be used to employ the structure of the differential equation under consideration to construct effective and robust methods. The global solution of PUFEM has been the theoretical basis of the local partition of unity finite element method, to be called later the extended finite element method.The first effort for developing the extended finite element method can be traced back to 1999, [7] presented a minimal enmeshing finite element method for crack growth. They added discontinuous enrichment functions to the finite element approximation to account for the presence of the crack. The method allowed the crack to be arbitrarily aligned within the mesh, though it required enmeshing for severely curved cracks. The method was improved and called extend Finite Element Method (XFEM) in [8, 9, 11, 14–16].In mathematics, relevant theoretical parts of doing are almost same in [17–19]. Of course, XFEM is still evolving currently and keeping the new extension about (such as fluid-structure interaction): mathematics stiff matrix remained prone to the sick.Consider a domainΩ in ℝ
n which is discredited by a set of elements Ω
e such that
(1)
⋃
e
∈
ε
Ω
e
=
Ω
,
and a set of nodes
(2)
𝒯
=
⋃
I
∈
𝒩
‍
x
I
.The total ordering of the element basis is given byε and the total ordering of the nodes is given by 𝒩. The standard finite element basis is defined as follows:
(3)
ℬ
std
=
⋃
I
∈
𝒩
‍
N
I
(
x
)
,
where N
I
(
x
) is the finite element basis or shape function of node I; N
I
(
x
) is assumed to be of compact support and piecewise continuously differentiable. Typically, ℬ
std spans the space of piecewise continuous polynomials of a specific order. The XFEM aims to alleviate the burden associated with mesh generation for problems with voids and interfaces.. It does not require the finite element mesh to conform to internal boundaries. The essence of the XFEM lies in subdividing a model into two distinct parts: mesh generation for the domain (excluding internal boundaries) and enrichment of the finite-element approximation by additional functions. A partitioning of a typical domain into its enriched and unenriched subdomains is shown in Figure 2(a). We defined the subdomain where the enrichment is applied as Ω
Enr This region’s feature is that the enrichment is dominant. The elements covering Ω
Enr are the enriched elements, so
(4)
⋃
e
∈
ε
Enr
‍
Ω
e
=
Ω
Enr
,
where ε
Enr is the ordering of the subset of elements which are to be enriched. Any node I for which the support of N
I
(
x
) overlaps or is contained in Ω
Enr is enriched node, and the set is defined by the ordering 𝒩
Enr, 𝒩
Enr
⊂
𝒩. For an enrichment function ψ
(
x
), the enriched basis for the local partition of unity method is
(5)
ℬ
Enr
=
ℬ
std
⊕
⋃
I
∈
𝒩
Enr
‍
N
I
(
x
)
ψ
(
x
)
.
The approximation enriched with a local partition of unity the field u
(
x
) is given by the sum of a standard finite element approximation and a linear combination of the enriched basis functions
(6)
u
h
(
x
)
=
∑
I
∈
𝒩
N
I
std
(
x
)
u
I
+
∑
J
∈
𝒩
Enr
N
J
Enr
(
x
)
ψ
(
x
)
q
J
,
where N
I
std
(
x
) and N
J
Enr
(
x
) are the shape functions for the standard part and the partition of unity, respectively; different order interpolates may be used for the standard and enriched shape functions.For example, we consider the linear elastostatics governing equations(7)
∇
·
σ
+
b
=
0
in
Ω
,
σ
=
C
:
ε
,
ε
=
∇
s
u
,
where σ is cauchy stress and b is the body force per unit volume. The first function is called equilibrium equation. The second function is the equilibrium equation that describes the relation of strain and displacement. The third equation is for Hooke’s law. The parameter C implies two parameters μ
=
G
=
E
/
(
2
(
1
+
ν
)
) and shear elasticity, λ
=
E
ν
/
(
1
+
ν
)
(
1
-
2
ν
) (Lamé constants are λ, μ and Poisson’s ratio is ν). Consider
(8)
u
=
u
¯
on
Γ
u
,
σ
·
n
=
t
on
Γ
t
,
σ
·
n
h
i
=
0
on
Γ
h
i
(
i
=
1,2
,
…
,
m
)
,
[
σ
·
n
I
i
]
=
0
on
Γ
I
i
(
i
=
1,2
,
…
,
n
)
.
The boundary Γ is composed of the sets Γ
u, Γ
t, Γ
h
i, and Γ
I
i such that
(9)
Γ
=
Γ
u
∪
Γ
t
∪
Γ
h
i
∪
Γ
I
i
,
where n is the unit outward normal to Ω and u
¯ and t are prescribed displacements and traction, respectively. The weak formulation of (12) is well known and reads as follows: find u
∈
V such that
(10)
∫
Ω
σ
(
u
)
:
ε
(
v
)
d
Ω
=
∫
Ω
b
·
v
d
Ω
+
∫
Γ
t
t
¯
·
v
d
Γ
∀
v
∈
V
0
.
Here V
=
H
1
(
Ω
) and V
0
=
{
v
∈
V
:
v
=
0
on
Γ
u
}.For the finite element discretization, letτ
h be a regular rectangular element of domain Ω and define the mesh parameter h
=
max
T
∈
τ
h
{
diam
(
T
)
}.Consider the Bubnov-Galerkin implementation for the XFEM in two-dimensional linear elasticity. In the XFEM, finite-dimensional subspacesV
h
⊂
V and V
0
h
⊂
V
0 are used as the approximating trial and test spaces. The finite dimensional space V
h is given by
(11)
V
h
=
{
v
∈
V
:
v
|
T
∈
Q
k
(
T
)
,
∀
T
∈
τ
h
}
.
The weak form of the discrete problem can be stated as follows: find u
h
∈
V
h such that
(12)
∫
Ω
h
σ
(
u
h
)
:
ε
(
v
h
)
d
Ω
=
∫
Ω
h
b
·
v
h
d
Ω
+
∫
Γ
t
h
t
¯
·
v
h
d
Γ
∀
v
h
∈
V
0
h
⊂
V
0
,
where V
h
⊂
V and V
0
h
⊂
V
0. In a Bubnov-Galkerkin procedure, the trial function u
h and the test function v
h are represented as linear combinations of the same shape functions. The trial and test functions are
(13)
u
h
(
x
)
=
∑
I
ϕ
I
(
x
)
u
I
+
∑
J
ϕ
J
ψ
(
φ
(
x
)
)
a
J
,
v
h
(
x
)
=
∑
I
ϕ
I
(
x
)
v
I
+
∑
J
ϕ
J
ψ
(
φ
(
x
)
)
b
J
,
where I denotes the set of all nodes in the mesh and J
=
{
j
∈
I
:
O
ε
x
j
∪
Γ
≠
∅
}
(
ε
≤
h
) denotes the set of nodes near the interface. ϕ
I
(
x
) are the finite-element shape functions, φ
(
x
) is the level set function, and ψ
(
φ
(
x
)
) is the enrichment function for material interface.There are several kinds of enrichment functions, such as abs-enrichment function, step-enrichment function, and Ramp-enrichment function. The choice of function is based on the behavior of the solution near the interface.The Ramp function is defined as(14)
ψ
Ramp
=
{
1
φ
>
0
1
-
2
φ
φ
≤
0
.
This enrichment function yields only continuous solutions. The advantage is that it automatically satisfies the continuity condition [
u
]
=
0 and does not require the use of Lagrange multipliers.The step-enrichment function is defined as(15)
ψ
Step
=
{
1
φ
>
0
,
-
1
φ
≤
0
.
This enrichment function can yield a continuous or discontinuous solution across the interface but requires Lagrange multipliers to apply the Dirichlet jump condition.On substituting the trial and test functions form (17), and using the arbitrariness of nodal variations, the following discrete system of linear equations is obtained:
(16)
K
d
=
F
,
where
(17)
K
I
J
=
∫
Ω
h
B
I
T
C
B
J
d
Ω
,
(18)
F
I
=
∫
Γ
t
h
ϕ
^
I
t
¯
d
Γ
+
∫
Ω
h
ϕ
^
I
b
d
Ω
.
ϕ
^
I
≡
ϕ
I for a finite-element displacement degree of freedom, and ϕ
^
I
≡
ϕ
I
ψ for an enriched degree of freedom. In the above equations, C is the constitutive matrix for an isotropic linear elastic material:
(19)
C
=
(
λ
+
2
G
1
0
1
λ
+
2
G
0
0
0
G
)
,
and the matrix B
I is given by
(20)
B
I
=
(
ϕ
^
I
,
x
0
0
ϕ
^
I
,
y
ϕ
^
I
,
y
ϕ
^
I
,
x
)
.We can get the integral terms, and the matrices computed from the integral terms are block matrices of the form(21)
K
=
(
K
U
U
K
U
A
K
A
U
K
A
A
)
,
F
=
(
F
U
F
A
)
,
d
=
(
U
Enr
)
.The next step is solving the linear equations system.


## 3. The Program for Solving the Equation

We need to solve the following equations:(22)
(
M
n
×
n
11
M
n
×
m
12
M
m
×
n
21
M
m
×
m
22
)
(
U
n
×
1
En
r
m
×
1
)
=
(
R
n
×
1
1
R
m
×
1
2
)
.
Here n
~
O
(
N
2
) and m
~
O
(
N
), the block matrix M
11 is same in FEM and (
M
21
)
′
=
M
12 are sparse matrix.We must point out that the matrixM
12 has up to m
2 nonzero elements. We can use the modified program, in (22) can be converted to two sub-problems: to solve (
M
22
)
-
1 and solve U by n-order linear equations. Now, we list the detailed steps.Algorithm 1.
(i) We use the ordinary method to solve(23)
M
11
U
(
0
)
=
R
1
,
the results for the FEM equation and set the result as U
(
0
), E
(
0
)
=
0. The matrix M
11 is block tridiagonal symmetric, sparse matrix. There are a lot of methods to solve this equation efficiently.
(ii) Solve the small scale equation(24)
M
22
En
r
(
i
+
1
)
=
R
2
-
M
21
U
(
i
)
,
i
=
0,1
,
…
.
The matrixM
22 is a sparse tridiagonal, where the dimension is much smaller than M
11. This equation can be solved quicker than (i).
(iii) Then we use (25) to solve U
(
i
+
1
);
(25)
M
11
U
(
i
+
1
)
=
R
1
-
M
12
En
r
(
i
+
1
)
,
i
=
0,1
,
…
.
,
the computation time complexity is the same with the first one.
(iv) If the norm∥
U
(
i
+
1
)
-
U
(
i
)
∥
≤
w or i
>
I
max (giving the maximum iterative times), output the result U
=
U
(
i
+
1
), Enr
=
En
r
(
i
+
1
). Else, return to step (ii).
In this part, we can use the common FEM solver such as LU or other useful solvers. We do not change the stiffness matrixM
11 and M
22, these two matrixes are symmetric, positive, definite and banded. We can compute easily for (23), (24), (25). The matrix conditions is not changed; 𝔎
(
M
11
)
=
O
(
h
-
2
) [17]. First of all, the store can choose a one-dimensional. We can choose a one dimensional storage when we use triangular decomposition for the matrix M.Remark 2.
If  the matrixM
11 and A
22 are not tridiagonal, we often use LU or GMRES method [20] to solve (23), (24), and (25).
The iteration has the following advantages.(1)
Other than the half bandwidth of zero elements do not have storage, nor involved in the calculation, the amount of storage crunch, the impact of rounding error is relatively small, even if it is ill-conditioned matrix, we often obtain a higher accuracy.
(2)
The operator only uses FEM linear equation solver; it is a simple calculation.
(3)
The original matrix deposit triangular matrix and diagonal matrix deposit area can overlap, thus avoiding to take up a lot of the middle unit of work, saving memory.
(4)
This solution can be extended to block direct solution.
(5)
With the iterative method, it is easy to build and test the general program and you can more accurately estimate the computation time; the user is relatively easy to grasp the algorithm process.Also if we using in [21] to solve the equation, (
(
2
N
)
2
+
4
N
) equations use LU solver which need O
(
(
2
N
)
2
+
4
N
)
3 Multiplication, here if the number of iterations is I, we use up to O
(
(
2
N
)
6
+
(
4
N
)
3
) multiplications.


## 4. Numerical Test

In this section, we report two experiments to verify the efficacy and accuracy of the new iterative method in XFEM. We use Matlab and modify software [22] to demonstrate Algorithm 1. For the  convenience of presentation, we introduce the following notations:SFEM means the standard finite element method;
XFEM means the extended finite element method;
h means the mesh size;
GMRES means the generalized minimal residual iterative method;
L
2-error means ∥
u
h
-
u
∥
0;
E
h means the L
2-error when the mesh size equals h;
iterations mean the iterative times in the algorithm;
orders mean the numerical error order by XFEM:(26)
orders
=
log
(
E
h
1
/
E
h
2
)
log
(
h
1
/
h
2
)
.Hereh
1
≤
h
2. First, we list Example 1 for one-dimensional problems.Example 1.
The exact solution is(27)
u
=
{
-
(
x
-
0.5
)
2
2
k
2
-
(
x
-
0.5
)
(
k
1
-
k
2
)
4
k
2
(
k
2
+
k
1
)
+
1
4
(
k
2
+
k
1
)
,
FDFDFDFDFDFDFDFDFDFDFDFiii
if
x
∈
Ω
1
,
-
(
x
-
0.5
)
2
2
k
1
-
(
x
-
0.5
)
(
k
1
-
k
2
)
4
k
1
(
k
2
+
k
1
)
+
1
4
(
k
2
+
k
1
)
,
ggdgdgdgdgdgdgdgdgdgdgdgdgdgdgg
if
x
∈
Ω
2
.HereΩ
1
=
[
0,0.5
]and Ω
2
=
(
0.5,1
] the coefficient k
=
1, ifx
∈
Ω
1, k
=
10, ifx
∈
Ω
2. And F is point force from the righ; F
=
10. We use Figure 1 to describe the L
2 error change when we use different h. We can easily see that, when the step is smaller than 1
/
80, the whole computing time has economized one order of magnitude in Example 1.(a) The comparison of time by iteration and LU with different mesh. (b) The comparison ofL
2-error by iteration and LU with different mesh.
(a)
(b)(a) All nodes and enrichment nodes (in red color);N
=
20; (b) distribution map of nonzero elements in whole matrix M; N
=
20.
(a)
(b)Example 2.
For the two-dimensional solid test, we use the same enrichment function and shape function to compute (4). We define Ω
=
[
0,1
]
×
[
0,1
] with the material interface
(28)
Γ
:
x
2
+
y
2
=
a
2
=
0.1608
.
Let Young’s modulus and Poisson’s ratio in
(29)
Ω
1
=
{
(
x
,
y
)
∈
Ω
∣
x
2
+
y
2
≤
0.404
2
}
be E
1
=
100, ν
1
=
0 and let those in
(30)
Ω
2
=
{
(
x
,
y
)
∈
Ω
∣
x
2
+
y
2
≤
1
4
}
be E
2
=
10, ν
2
=
0, but area force in domain f
x
=
f
y
=
0.
We denoter
=
x
2
+
y
2 and tan
θ
=
y
/
x, and the exact solution is
(31)
u
=
a
8
μ
(
-
2
a
3
r
3
cos
3
θ
)
r
a
(
κ
+
1
)
cos
θ
+
2
a
r
(
(
1
+
κ
)
cos
θ
+
cos
3
θ
)
g
f
g
f
g
f
g
g
-
2
a
3
r
3
cos
3
θ
)
,
v
=
a
8
μ
(
-
2
a
3
r
3
sin
3
θ
r
a
(
κ
-
3
)
sin
θ
+
2
a
r
(
(
1
-
κ
)
sin
θ
+
sin
3
θ
)
g
f
g
f
g
f
g
-
2
a
3
r
3
sin
3
θ
)
.Plate with a circular hole. Figure2(b) shows the consistent rectangular subdivision and labels the nodes which need to be extended. The CPU time we counted is the total time for the entire programm running. Table 1 shows the linear equation scale and the condition numbers of FEM and XFEM, the density of sparse matrix show the matrix is sparse. It can be easily find the condition number of XFEM is bigger than that of FEM. The CPU times are also compared with LU solvers and GMRES in Table 2. The cheapest is LU by FEM and the new iterative method is the second cheapest; however, the error for LU by FEM is the biggest in all methods. Table 3 shows an L
2-error comparison between the new iterative method and other solvers. And the numerical L
2-error is nearly O
(
h
2
). L
2-error, order, and iterations are given in Table 4.Table 1
DOF, sparsity of matrix, and condition number of linear equation by FEM and XFEM for Example2.
h
DOFFEM
DOFXFEM
condsFEM
condsXFEM
Sparsity of matrix (%)
1
20
882
1018
1.7627
×
1
0
4
7.90617
×
1
0
4
1.72
1
40
3362
3626
3.75035
×
1
0
4
3.20917
×
1
0
5
0.4774
1
60
12800
13120
1.02381
×
1
0
5
1.30023
×
1
0
6
0.222
1
80
12800
13120
1.02381
×
1
0
5
1.30023
×
1
0
6
0.13
1
100
25600
26240
4.07606
×
1
0
5
5.23974
×
1
0
6
0.0829Table 2
The CPU time for different solvers in Example2.
h
FEMLU
XFEMiteration
XFEMLU (s)
XFEMGMRES (s)
1
20
1.2166
1.3129
2.6019
2.1449
1
40
5.7133
6.5198
7.3678
8.9494
1
60
17.8348
20.6117
27.2810
26.6332
1
80
54.3007
59.0695
70.4340
70.7011
1
100
137.89
145.92
165.1256
166.5808Table 3
The comparison ofL
2-error with different solvers for Example 2.
h
FEMLU
XFEMLU (s)
XFEMGMRES
XFEMiteration (s)
1
20
4.2101
×
1
0
-
2
6.4292
×
1
0
-
3
6.3756
×
1
0
-
3
6.3756
×
1
0
-
3
1
40
2.1032
×
1
0
-
2
1.7108
×
1
0
-
3
1.7097
×
1
0
-
3
1.7092
×
1
0
-
3
1
60
1.4089
×
1
0
-
2
7.7599
×
1
0
-
4
7.7595
×
1
0
-
4
7.7554
×
1
0
-
4
1
80
1.0317
×
1
0
-
2
4.1781
×
1
0
-
4
4.1434
×
1
0
-
4
4.1733
×
1
0
-
4
1
100
8.2802
×
1
0
-
3
2.7344
×
1
0
-
4
2.7002
×
1
0
-
4
2.7295
×
1
0
-
4Table 4
L
2-error, order, CPU time, and iterations by XFEM with the new iterative method for Example 2.
h
ErrorXFEM
Orders
CPUXFEM (s)
Iterations
1
20
6.3754
×
1
0
-
3
∖
1.313
24
1
40
1.7092
×
1
0
-
3
1.8992
6.52
19
1
60
7.7554
×
1
0
-
4
1.9489
20.61
19
1
80
4.1733
×
1
0
-
4
2.1541
59.07
19
1
100
2.7295
×
1
0
-
4
1.9028
145.92
18


## 5. Conclusion

In this paper, we discuss how to use XFEM for the stiff matrix ill-conditioned and came to the following conclusions. Firstly, we found why some of problems cannot be computed by XFEM; the stiff matrix is not then symmetric and ill-conditioned. So, a lot of storage space is wasted. The time of computing XFEM equation is longer, which cannot be computed in certain subdivision. Secondly, we use XFEM for solving two-dimensional solid test; when solving the XFEM equations, the iterative method and the standard finite element equations solver are adopted, so that we change the problems into solving the smaller scale equations and standard finite element stiffness matrix need not be changed. Finally, numerical examples show that the proposed method is more efficient than the common method. We can use this programm in many other interface problems later. Besides, we also consider combining the stable GFEM with our iterative method in the future. And it is believed that the XFEM can be widely used in more problems in fluid dynamics and material science. The calculation module in some software can be improved for users.


## Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.



---

*Source: 367802-2014-01-14.xml*