XLME
interpolants,
a seamless
bridge
between
XFEM
and
enriched
meshless
methods
F.
Amiri
• C.
Anitescu
•
M.
Arroyo
• S.
P.
A.
Bordas
•
T.
Rabczuk
Abstract
In
this paper.
we
develop a method based
on
local
maximum entropy shape functions together with enrichment
functions used
in
partition
of
unity
methcx:ls
to discretize
problems in linear elastic fracture mechanics. We obtain
improved accuracy relative to the standard extended finite
element method
at
a comparable computational cost.
In
addi-
tion.
we
keep the advantages
of
the LME shape functions.
such as smootlmess and non-negativity. We show numeri-
cally that optimal convergence (same as in FEM) for energy
norm and stress intensity factors can be obtained through
the
use
of
geometric (fixed area) errriclnnent
with
no
special
treatment of the nodes near the crack such as blending or
shifting.
Keywords Local maximum
entropy·
Convex
approximation· Meshless
methods·
Extrinsic enrichment
S.
P.
A. Bordas's ORCID ID is 0000-0001-7622-2193.
F.
Amiri . C. Anitescu .
T.
Rabczuk
(~)
Institute
of
Structural Mechanics, Bauhaus University Weimar,
Marienstr.
15,99423
Weimar,
Gennany
e-mail: timon.rabczuk@uni-weimar.de
F.
Amiri
e-mail: fatemeh.amiri@uni-weimar.de
M. Arroyo
Departament de
Matematica Aplicada
3,
School
of
Civil
Engineering
of
Barcelona (ETSECCPB), Universitat
Politecnica de Catalunya, Barcelona, Spain
S.
P.
A. Bordas
Cardiff School
of
Engineering, Institute
of
Mechanics
and Advanced Materials, Cardiff University, Queen's Buildings,
The
Parade, Cardiff, Wales CF24
3AA,
UK
T.
Rabczuk
School
of
Civil, Environmental and Architectural Engineering,
Korea University, Seoul, South Korea
1
Introduction
Maximum entropy shape functions are a relatively new class
of
approximation functions, as they were first introduced in
[1]
in
the context
of
polygonal interpolation.
The
idea
of
these
functions is to maximize the Shannon entropy [2]
of
the
basis
functions, which gives a measure of
the
uncertainty
in
the
approximation scheme.
The
principle
of
maximum entropy
(max-ent) was developed
by
Jaynes [
3.4
].
who
showed that
there is a natural correspondence between statistical mechan-
ics and information theory.
In
particular, max-ent offers the
least-biased statistical inference
when
the
shape functions
are viewed as probability distributions subject to the
approx-
imation constraints (such as linear reproducing properties).
However, without additional constraints,
the
basis functions
are non-local, which due to increased overlapping makes
them unsuitable for analysis using Galerkin methods.
The
increased overlapping
of
the
basis functions generally leads
to more expensive numerical integration due to
the
largenum-
ber of evaluation points.
It
also produces a non-sparse stiff-
ness matrix, resulting
in
a linear system that is
much
more
expensive to solve.
The
local maximum-entropy (LME) approximation
schemes
were
developed
in
[5] using a framework similar
to meshfree methods. Here
the
support
of
the basis functions
is introduced as a thermalization (or penalty) parameter
fJ
in
the
constraint equations.
When
f3
= 0, then
the
max-
ent principle is fully satisfied and
the
basis functions will
be
least biased. For example,
if
only zero-order consistency
is required, the shape functions are Shepard approximants
[6]
with
Gaussian weight function.
When
fJ
is large. then
the shape functions have minimal support.
In
particular. they
become the usual linear finite element functions defined on
a Delaunay triangulation
of
the
domain associated
with
the
given node set.
In
[5] it was shown that for some values
of
fJ.
1
the approximation properties of the maximum-entropy basis
functions are greatly superior to those
of
the finite element
linear functions, even when the added computational cost
due to larger support is taken into account.
Subsequent studies, such as [7- 9
],
show that maximum
entropy shape functions are suitable for solving a variety
of
problems such as thin shell analysis, compressible and
nearly-incompressible elasticity and incompressible media
problems. Higher order approximations can also be obtained
using the max-ent framework, as shown in [10
].
This class
of
methods is therefore related to the MLS-based meshless
methods (due to the node-based formulation) and isogeo-
metric analysis (with whom it shares features such as weak
Kronecker delta and non-negativity), inheriting some advan-
tages from both.
In this work, we propose a coupling
of
the LME shape
functions with the extrinsic enrichments used in partition
of
unity enriched methods for fracture, such as the extended
finite element method (XFEM), see [
11
-
13
].
There
is
a growing interest in modeling fracture mechan-
ics with enrichment functions combined with meshless meth-
ods [14- 16
],
isogeometric analysis [
17
],
or strain-smoothed
XFEM [18,
19
].
Advantages
of
the meshless and isogeomet-
ric methods include the possibility to model curved bound-
aries through higher order shape functions and to resolve the
gradient fields more accurately than with low order Lagrange
elements. This higher regularity
of
the basis functions
is
also
particularly advantageous when the model problem requires
it,
such as for the Kirchhoff-Love theory. Also in some
enriched meshless methods, no representation
of
the crack's
topology is needed as this is handled through cracking parti-
cles as in [20] or weight-function enrichments as in [
21
,
22
].
Here, we show that the enriched maximum entropy shape
functions are suitable for this class
of
problems. Moreover,
this method
is
more accurate than standard XFEM and does
not require the so-called blending elements (the elements
near the crack tip). When compared
to
usual meshfree meth-
ods for crack propagation, such as Element Free Galerkin
(EFG), the method presented here can more easily deal with
essential boundary conditions, due to the fact that the shape
functions satisfy a weak Kronecker delta property. The shape
functions are also very smooth
(COO),
which results in an
accurate numerical integration with a relatively low number
of
integration points, especially for Gauss-Legendre quadra-
ture [5, 8, 10
].
Moreover, smooth and non-negative basis func-
tions, such as those used in isogeometric analysis are gaining
impetus.
The paper is organized as follows: in the next section we
will briefly describe the LME approximants. Then we will
introduce the coupling between LME and XFEM, with par-
ticular reference to implementation issues such as numerical
integration. Next we examine the accuracy
of
the method
through several numerical examples, which indicate that the
convergence rates for the energy norm
of
the error and the
stress-intensity factors, are
O(h)
and O(h2) respectively.
Some concluding remarks are stated in the last section.
2 Local
maximum
entropy
(LME)
approximants
Local maximum entropy meshfree approximants, introduced
in [5], are related to other convex approximation schemes,
such as natural neighbor approximants [23
],
subdivision
approximants [24
],
or B-spline and NURBS basis func-
tions [25
].
The LME basis functions will be denoted
by
Pa
(x), a =
1,
...
, N with x E
lR:
d
,
d is the dimension
of
the physical domain. They are non-negative and are required
to satisfy the zeroth-order and first-order consistency condi-
tions:
Pa
(x) 2' 0,
(1
)
N
LPa(X)
=
1,
(2)
a=l
N
L
Pa
(x)xa = x.
(3)
a=l
In
the last equation, the vector
Xa
identifies the positions
of
the nodes associated with each basis function. Consider a set
of
nodes X =
{Xala~I,
...
,N,
which we will call the node set.
The convex hull of X is the set
convX
:=
{x
E
lR:dlx
=
XA,
A E
lR:~,
1·
A =
1)
(4)
Here]R~
is the non-negative orthant, 1 denotes the vector in
]RN
whose entries are one, and X is the d x N matrix whose
columns are the co-ordinates
of
the position vectors
of
the
nodes in the node set X [5
].
Convex approximants, which are
in the span
of
convex basis functions, can only exist within
the convex hull
of
X (or subsets
of
it) and satisfy a weak Kro-
necker delta property at the boundary
of
the convex hull
of
the nodes. This means that the shape functions correspond-
ing to the interior nodes vanish on the boundary. With this
property, the imposition of essential boundary conditions in
the Galerkin method is straightforward.
The principle of maximum entropy comes from statisti-
cal physics and information theory, which consider the mea-
sure of uncertainty or information entropy [2]. Consider a
random variable
X I
--+
]Rd, where I
is
the index set
I = {I,
...
,
N)
andx(a)
=
Xa
gives to each index theposi-
tion vector
of
its corresponding node. Since the shape func-
tions
of
a convex approximation scheme are non-negative
and add to one, we regard
{PI
(x),
...
, PN(X)) as the corre-
sponding probabilities. The statistical expectation or average
of
this random variable, as regarding Eq. (3
),
is x. Accord-
ing to this interpretation, the approximation
of
a function
u(x)
'"
L~~IPa(X)Ua
from the nodal values
{Uala~I,
...
,N
2
is
understood as an expected value
u(x)
of
a random variable
!"
:
1-+
lR:
where !"(a) = U
a
.
The main idea
of
max-ent is to maximize the Shannon's
entropy,
H (PI,
P2,
...
, PN), subjectto the consistency con-
straints as follows:
(ME)
For a fixed x maximize
N
H(PI,
P2,···,
PN) = -
LPa
10g(Pa)
a=l
subject to
Pa
C>
0, a =
1,
... , N
N
LPa
= 1
a=l
N
LPaXa
= X
a=l
(5)
Solving the
(ME)
problem produces the set
of
basis functions,
Pa
:=
Pa(X), a =
1,
...
,
N.
However, these basis functions
are non-local, i.e. they have support in all
of
conv X, and are
not suitable for use in a Galerkin approximation because it
would lead to a full, non-bauded matrix. Nevertheless, they
have been used in
[1] as basis functions for polygonal ele-
ments.
Another optimization problem which takes into account
the locality of the shape functions is Rajau's form of the
Delaunay triangulation [26]. This can be stated as the fol-
lowing linear program:
(RAJ) For a fixed x minimize
N
U(X,PI,P2,···,PN)
=
LPaI
X
-
X
aI
2
a=l
subject to
Pa
C>
0, a =
1,
...
, N
N
LPa
= 1
a=l
N
LPaXa
= X
a=l
(6)
ItiseasytoseethatU(x,
PI,
P2,
...
,
PN)isminimizedwhen
the shape functions PI,
...
,
PN
decay rapidly as the distauce
from the corresponding nodes
Xa
increases. There, the shape
functions that satisfy (RAJ) problem will have small sup-
ports, where the support can be defined up to a small tolerauce
E by
sUPP(Pa)
= {x :
Pa
(x) >
E)
The main idea
of
LME approximants
is
to compromise
between the
(ME) problem aud the (RAJ) problem
by
intro-
ducing parameters
fJa
that control the support
of
the
Pa.
Therefore we write:
For a fixed x minimize
(7)
N N
LfJaPa
Ix
- xal
2
+
LPa
10g(Pa)
a=l a=l
subjectto
Pa
C>
0, a =
1,
...
, N
N
LPa
= 1
a=l
N
LPaXa
= X
a=l
The non-negative parameters
f3a
can in general be functions
of
the position x. This convex optimization problem is solved
efficiently by a duality method as described in
[5]. Finally,
the shape functions are written in the form:
1
Pa(X)
=
exp[-fJa
Ix-x
a
I
2
+A*(X)'
(x - xa)]
Z(x,
A*(X»
where
N
Z(X,A) =
Lexp[-fJblx-xbI2+A'
(X-Xb)]
b=l
is
a function associated with the node set X and A *(x) is
defined by
A*(X)
= arg min log
Z(x,
A)
AElRd
The local max-ent shape functions are as smooth as fJ(x)
aud
Pa
(x,
fJa)
is a continuous function
of
fJ
E [0,
+00)
[5
].
For example LME shape functions are
Coo
if
fJ
is constant.
In
this paper we choose
f3
=
fz,
where h is a measure
of
the
nodal spacing and y is constant over the domain.
In
this case
the shape functions are smooth and their degree
of
locality
is
controlled by the parameter
y.
A plot
of
the LME functions
for y
= 1.8 and a particular choice
of
nodes is given in
Fig.
!.
In general, the optimal
fJ
is not obvious aud this will
be discussed later in this paper.
As we mentioned before, LME shape functions satisfy a
weak Kronecker delta property at the boundary of the con-
vex hull
of
the nodes. Therefore, the shape functions that
correspond
to
interior nodes vanish on the boundary.
3 Brief on extrinsic enrichments for partition
of
unity
methods
3.1
Description
The main idea of partition
of
unity (PU) enrichment as used
here is to extend the max-ent approximation space with some
additional enrichment functions. The proposed method is
based on a local
PU
and uses an extrinsic enrichment to
3
0.8
0.6
0.4
0.2
O.
•
0
2
•
• •
8 8
6
•
•
4
Fig. 1 Local max-ent shape functions in 2D
•
•
2
model the discontinuity. The max-ent approximation can
be
decomposed into a standard part and an enriched part:
uh(x) = L PI (X)UI + L
PI(x)X(1o(x»aJ
lEW
4
+ L PK(X) L Bk(X)bkK
KEWs
k=l
Here the first term is the standard approximation part and
the second and the third terms are the enriched parts. W
is
the set
of
nodes in the entire discretization and
Wb
and
Ws
are the sets
of
enriched nodes. P I are the shape func-
tions and X and
Bk
are the enrichment functions. Normally,
X is selected as a step or Heaviside function and is used
to enrich the nodes where the supports
of
the
LME
shape
functions are completely cut by the crack.
Bk
are branch
functions and are used to errrich the shape functions whose
supports include the crack tip.
In
this paper
we
use a geomet-
ric (fixed area) enrichment. and therefore
we
obtain optimal
convergence rate
[0
(h
2
)]
without a special treatment
of
the
so-called
"blending" area around the crack tip. Branch func-
tions are defined as follows
(in polar coordinate relative to
the crack tip. denoted by
x
tip
):
. e
Bl
(r, e) =
~
sm
2:
e
B2(r, e) =
~
cos
2:
. e
B3(r, e) =
~
sm
2:
cose
e
B4(r, e) =
~
cos
2:
cose,
where r =
Ix
- xtipl·
(8)
(9)
(10)
(11)
10
(x)
is
the signed distance from the point x to the crack
segment and
aI
and
bkI
are additional degrees
of
freedom
[27].
The
signed distance function
is
defined as:
x
r
n
¢<
o
¢>o
Fig. 2 Signed distance function
1o(x)
= min Ix -
xrl
sign(n·
(x -
xr»
xr
Er
Here r is the curve
of
discontinuity,
Xr
is
an
arbitrary point
on
rand
n is normal vector to r (see Fig. 2
).
If
we
choose
X as a Heaviside function, then
H(1o(x» =
{I
-1
if
10
(x) > a
if
10
(x) < a
(12)
This enrichment function captures the
jump
across the crack
faces.
In
order to model a curved crack, the signed distance
function can
be
approximated by the same shape functions
as the displacement. Assume
t is a vector tangent to the
curved crack, directed towards the crack tip. We approximate
10
by:
(13)
Here
10
I are the nodal values
of
10,
PI are the shape functions
and
Q¢,
is the domain
of
definition for
10,
given by:
Q¢
:=
{xlt·
vr(x)
> 0) (14)
So, the approximated crack position
is
considered as:
r:=
{xl¢(x)
=
o,x
E
Q¢)
(15)
In
this case,
¢(x)
is not defined beyond the crack tip. So,
&0
possibilities are considered for the angle e
of
the Branch
functions.
If
t·
v r S
0,
then the regular polar angle from
-t
is
computed.
If
t .
vr
> 0, e is considered as in [28
]:
-10
e = arctan( )
Jr2
_
10
2
(16)
4
3.2 Numerical integration
3.2.1 Numerical integration for
IME
The
numerical integration
of
LME shape functions poses
similar challenges as that
of
the shape functions used in mesh-
less methods.
In
particular, the integrands used in the assem-
bly
of
the stiffness matrix are non-polynomial and (depend-
ing
on
the values
of
the parameter
y)
the supports of the shape
functions overlap more than in standard finite elements. How-
ever, the shape functions are smooth so only a relatively small
number
of
integration points are required.
In
the examples
we
considered,
we
used quadrilateral
background integration cells for integrating the shape func-
tions whose support does
not
intersect the crack. For the
values
of
y between 4.8 and 1.8, and for uniformly spaced
nodes and square
we
fouud that the 4 x 4 Gauss quadrature
rule is sufficient
to
ensure optimal convergence. Moreover,
a quadrature rule with
8 x 8 Gauss points provides close to
exact integration (i.e. the results change by less than
10-
6
when
the number
of
Gauss points is further increased).
3.3 Numerical integration for enriched LME
The
usual numerical integration methods, for example Gauss
quadrature, are less accurate for
PU-enriched methods for
fracture. This happens due to the discontinuity along the
crack,
and
the singularity at the crack tip.
The
usual rule
is to use a simple splitting
of
integration cells crossed by the
crack [29
].
In
[30], a method was proposed in which each
part
of
the elements that are cut or intersected by a discontinuity
is mapped onto the unit disk using a conformal
Schwarz-
Christoffel map. However, for straight cracks, a triangulation
of
the elements
cut
by the crack which takes into account the
location
of
the discontinuity is relatively easy to implement
and was used in this work.
For the integration cells that contain the crack tip, special
care has to be taken. These cells contain the discontinuity
and
a singularity together. So, simply refining the triangles that
make up the integration cells leads to less accurate numerical
results. A simple solution is to refine locally each split trian-
gle, until
an
acceptable estimate
of
the integrands is achieved.
Unfortunately, this method is expensive.
To
solve this prob-
lem, the
almost polar integration was introduced in [
29
].
The
main idea is to build a quadrature rule
on
a triangle from
a quadrature rule
on
the unit square (see Fig. 3
).
The
map
is:
T :
(x,
y)
--+
(xy,
y)
which maps a square into a triangle. By looking at the inte-
grands which contain the derivatives
of
the branch functions,
we
notice that the Jacobian
of
the transformation T, will
cancel the
r-
1
/
2
singularity. This integration method gives
_
crac
k tip
o
•
1
2
12
..
3
•
4
o
3
•
Fig. 3 Transfonnation of an integration method on a square into
an
integration method
on
a triangle for crack tip functions
excellent results
with
a low number of integration points and
is used
on
the sub-triangles having the crack tip as a vertex.
In
the other integration cells,
we
found it is sufficient to use
standard Gauss quadrature over a background mesh (such
as the Delaunay triangulation
of
the nodes that takes in to
account the discontinuity for the cells cut by the crack).
An
important distinction between meshless methods and
standard finite elements is that, in the former, the numeri-
cal integration is almost never exact. Recent work
[
31
] has
shown that integration errors in meshless methods negatively
impact the stability of the method
when
a large number
of
degrees of freedom is involved.
In
particular, as the value
of
the discretization parameter h decreases, the accuracy
of
the
numerical integration should increase proportionally, so that
optimal convergence can be obtained. We have conducted a
ueLaileu
sLuuy
on
Ule
dfed
of
approxilllaLe inLegraLion for
one of the numerical examples shown below.
3.4 Condition number
There are two ways to choose the enrichment area:
topologi-
cal enrichment
in which the area
of
enrichment shrinks with
the nodal spacing
h,
and geometric enrichment which uses a
fixed enrichment area.
In
topological enrichment, the branch
functions are multiplied by shape functions
on
a small set
of
nodes around the crack tip. These singular functions live on
a compact support vanishing as h goes to zero.
In
the con-
text of meshless methods, only topological enrichment has
been studied, which leads to non-optimal convergence rate.
However, the numerical results
of
this paper show that the
enrichment area should have a size independent
of
the mesh
parameter (i.e. it should be geometric) to obtain optimal con-
vergence, as seen for standard XFEM in [29,32
].
Unfortu-
nately, adding singular functions on all the nodes within a
fixed area around the crack tip leads to an increase in the
number
of
degrees
of
freedom and an increase in the condi-
tion number (see Fig. 4
).
Some methods were proposed to improve the condition
number of the stiffness matrix, such as preconditioning
schemes. Here
we
use a method introduced in [32] which
relies
on
a Cholesky decomposition
of
the diagonal blocks
5