Adaptive finite element analysis of elliptic problems based on bubble-type local mesh generation

Accepted Manuscript Adaptive finite element analysis of elliptic problems based on bubble-type local mesh generation Weiwei Zhang, Yufeng Nie, Yuantong Gu PII: DOI: Reference:

S0377-0427(14)00528-7 http://dx.doi.org/10.1016/j.cam.2014.11.049 CAM 9892

To appear in:

Journal of Computational and Applied Mathematics

Received date: 1 April 2014 Revised date: 30 September 2014 Please cite this article as: W. Zhang, Y. Nie, Y. Gu, Adaptive finite element analysis of elliptic problems based on bubble-type local mesh generation, Journal of Computational and Applied Mathematics (2014), http://dx.doi.org/10.1016/j.cam.2014.11.049 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript Click here to view linked References

Adaptive finite element analysis of elliptic problems based on bubble-type local mesh generation Weiwei Zhanga , Yufeng Niea,∗, Yuantong Gub a Department

of Applied Mathematics, School of Science, Northwestern Polytechnical University, Xi’an 710129, PR China b School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology, Brisbane, Australia

Abstract A new mesh adaptivity algorithm that combines a posteriori error estimation with bubble-type local mesh generation (BLMG) strategy for elliptic differential equations is proposed. The size function used in the BLMG is defined on each vertex during the adaptive process based on the obtained error estimator. In order to avoid the excessive coarsening and refining in each iterative step, two factor thresholds are introduced in the size function. The advantages of the BLMG-based adaptive finite element method, compared with other known methods, are given as follows: the refining and coarsening are obtained fluently in the same framework; the local a posteriori error estimation is easy to implement through the adjacency list of the BLMG method; at all levels of refinement, the updated triangles remain very well shaped, even if the mesh size at any particular refinement level varies by several orders of magnitude. Several numerical examples with singularities for the elliptic problems, where the explicit error estimators are used, verify the efficiency of the algorithm. The analysis for the parameters introduced in the size function shows that the algorithm has good flexibility. Keywords: Mesh size function; node placement; local mesh generation; error estimator; adaptive analysis

∗ Corresponding

author. Tel.: +86 29 88431606; fax: +86 29 88431655. Email address: [email protected] ( Yufeng Nie)

Preprint submitted to Elsevier

September 30, 2014

1. Introduction For solving elliptic boundary value problems for which the solution has singularities, adaptive finite element methods (AFEMs) have the potential advantage of a significant reduction of the computational cost [1–3], through designing reliable and efficient a posteriori error estimators and mesh adaptive techniques. Both mesh adaptivity and a posteriori error estimators have been extensively studied, beginning with the pioneering work of Babuˇska etc. [4], and then developed by many others researchers [5–9] and those referenced therein. The performance of adaptive methods for PDEs depends not only on the error estimators, but also on the techniques used for adaptively generating meshes [10]. In high quality adaptive mesh generation, different methods exist to perform mesh adaptation and they can be broadly classified into three categories: h-type refinement, p-type refinement and r-type refinement. In the h-type refinement, the meshes are coarsened in regions where errors are relatively small and refined in regions where the errors are relatively large. Two families exist in this class of methods [11]. In the first one, a new mesh is built according to a required mesh size. For example, the centroidal voronoi Delaunay triangulations have been developed for adaptive mesh generation and optimization by Ju [10] and Huang [12]. In the second family, a new mesh is obtained by dividing some selected elements. One common subdivision algorithm for the selected elements is the bisection method [13, 14] (the longest side bisection method and the newest vertex bisection method). This approach is very simple to implement, however, the mesh smoothing is usually necessary for generating unstructured mesh, and some appropriate measures are needed to remove the hanging nodes, coordinating the new mesh. Another subdivision algorithm is to use regular grids with the 1-irregularity rule, which allows the presence of the hanging nodes, see Refs. [15, 16]. While some difficulties may exist for this subdivision algorithm to partition the domains with complex geometric boundaries. Meanwhile, an irregular high-order polynomial interpolation caused by the hanging nodes is required in this subdivision algorithm for the finite ele-

2

ment solution. In some sense it increases computational complexity compared to using linear interpolation polynomial. In the p-type refinement, the same mesh is used but the order of the polynomial functions is increased. Since that implies large modifications in the simulation software, and the improvement of accuracy when using high order polynomial functions is largely based on the smoothness of finite element solution, they are not widely used. In the r-type refinement, the connectivity of the mesh is kept but the nodes are moved to increase mesh density in some locations. Since new nodes are not added, the accuracy of the solution derived by this method is limited by the initial number of nodes and elements. Moreover, the three categories are often combined into hybrid refinement strategies, such as, hp-type [17–19] and hr-type [20] adaptive methods. In this paper, the focus is on the adaptive finite element analysis based on bubble-type local mesh generation (BLMG) method. To be precise, the BLMG-based adaptive method is an h-type unstructured refinement method, compared with the h-type adaptive method with regular grids, it has better adaptability to partition the domains with complex geometric boundaries, and could be used in the situation of anisotropic case. Initially, Shimada et al. [21] proposed the bubble packing methods (BPM) to determine the nodal positions, and Delaunay triangulation technique is used to generate mesh. Then through combining Shimada’s method and Zienkiewicz and Zhu error estimator [22], Kim et al. [23] developed an adaptive bubble packing technique which leads to a reduction of errors with the increase of the refinement level. Moreover, Chung et al. [24] developed a remeshing algorithm using the BPM and the recoverybased error estimation approach [25], and applied the algorithm to the metal forming problem. In order to avoid the situation that the movement, addition or deletion of nodes coexist with mesh topology, a pure node placement method by bubble simulation was proposed by Nie et al. [26]. Instead of utilizing Delaunay triangulation to define the neighborhood of a bubble topologically, an adjacency list is introduced [26] for each bubble which stores the information of neighbor 3

bubbles. Then the parallel version [27] and a series of acceleration strategies [28] for this node placement method have been developed subsequently. Further, based on the generated high-quality nodes set and the information of adjacency lists, the node-based local mesh generation method BLMG has also been presented in [29]. In the BLMG method, the series of acceleration strategies, the improvement of viscosity coefficient in motion equation of bubbles and the node connection algorithm make the BLMG completely different from the above-mentioned BPM. The BLMG method has been applied for anisotropic triangular meshing [30] and parametric surface triangulation problems [31]. While the use of BLMG in the AFEM context has not been previously explored, thus in this study, we aim to develop a convergent AFEM with the BLMG for two dimensional elliptic problems based on error estimation. A brief outline of our method is now given. After a coarse Delaunay triangular mesh with a uniform size is completed, the iterative loop is executed which consists of computing the finite element solution with respect to the current partition, computing an a posteriori error estimator, defining the new size function and the construction of the next partition. At each refinement level, based on the obtained error estimator, a new value of mesh size is calculated on each vertex and the corresponding mesh size field function is obtained by linear interpolation. Then the nodes are gradually moved into a near-optimal configuration by performing dynamic simulation based on Newton’s second law of motion. At last, the BLMG method is applied to form the patch of elements around each node, and the union of the patches is a Delaunay triangulation of the domain. The updated triangular mesh is of high quality so as to that the adaptive process can achieve high-accuracy finite element solutions with a small degree of freedom. The remainder of this paper is organized as follows. In Section 2, we review the specific a posteriori error estimators used in this work. In Section 3, how to modify the mesh size function is firstly shown according to the a posteriori error estimators, where the size function is used to decide how mesh nodes should be distributed. Then the node placement method by bubble simulation is briefly 4

described, and the details of the BLMG method are introduced later. In Section 4, the flow chart of our mesh adaptation algorithm is given, and some numerical examples are presented in Section 5 to check the validity and efficiency of the present approach. At last, conclusions and future works are given in Section 6. 2. Error estimator for linear finite element methods Assume that Ω ⊂ R2 be a bounded domain with a Lipschitz boundary ∂Ω.

T is a triangulation of Ω. We consider the following second order elliptic partial differential equation

  

−∆u = f in Ω

u = 0 on ∂Ω

(1)

where f is a smooth function defined in Ω. For effective and robust adaptation of the finite element solution procedure, it is necessary to derive a reliable error estimator in terms of the computed finite element solution uh . There are several types of a posteriori error estimators used in AFEMs [7, 8, 32], e.g., explicit error estimators, implicit error estimators and the recovery-based error estimators. Explicit error estimators [4] involve a direct computation of the interior element residuals and the jumps at the element boundaries to find an estimate for the error in the energy norm. Implicit error estimators [5, 9] require the solution of auxiliary local boundary value problems. The recovery-based error estimators [12, 13, 22] smooth the gradients of the solution and compare the unsmoothed and the smoothed gradients in order to assess the solution error. Here, the construction of a posteriori error estimators is not our focus. Two conventional explicit error estimators are used to detect local regions of large error, and they are combined with the adaptive mesh refinement method BLMG; the effectiveness and the convergence properties of our method are demonstrated later. Let V be the Hilbert space H01 (Ω). The weak form of the above model problem (1) is: find u ∈ V such that a(u, v) = l(v) ∀v ∈ V, 5

where the bilinear form a(u, v) and the linear functional l(v) are defined as Z Z a(u, v) = ∇u · ∇vdΩ, l(v) = f vdΩ. Ω

Ω

For any u ∈ V , we define its energy norm || · ||E by ||u||E =

p a(u, u).

Let Vh ⊂ V be the space of continuous piecewise linear finite element space associated with the triangulation T . The finite element approximation is: find

uh ∈ Vh such that

a(uh , vh ) = l(vh )

∀vh ∈ Vh .

Let eh = u − uh be the error of the approximate solution uh , satisfying the

error representation

a(eh , v) = a(u, v) − a(uh , v) = l(v) − a(uh , v)

∀v ∈ V.

(2)

If the domain integral is split into the contributions from each element, (2) can be rewritten as Z X Z a(eh , v) = ∇uh · ∇vdΩ), ( f vdΩ − T ∈T

(3)

T

T

where T denote an element in T . Let εI denote the set of interior edges of T , if

T and T ′ share the common edge γ ∈ εI , then the jump of the gradient across

the edge γ is

[(∇uh ) · nγ ] = (∇uh )|T · nT + (∇uh )|T ′ · nT ′ , where nT is the unit outward normal vector to ∂T . Applying integration by parts to the second term in (3) and rearranging terms, we can get XZ XZ a(eh , v) = RvdΩ + Jvds ∀v ∈ V, T ∈T

T

γ∈εI

(4)

γ

where R = f + ∆uh is the interior element residual, and J = −[(∇uh ) · nγ ]

defined on the edge γ represents the jump discontinuity in the approximation to the normal flux on the interface. 6

For given v ∈ V , let Ih v be the interpolant of v in Vh . Then, utilizing the

Galerkin orthogonality property a(eh , Ih v) = 0 and (4), we have a(eh , v) =

P R

T

T ∈T

R(v − Ih v)dΩ +

P R

γ∈εI

γ

J(v − Ih v)ds

∀v ∈ V

(5)

Applying the Cauchy-Schwarz inequality elementwise and the interpolation theory, and substituting eh in place of v, we can obtain the a posteriori error estimate

||eh ||2E

(

P

P

+ γ∈εI n o P 1 2 2 hT ||R||L2 (T ) + 2 hT ||J||2L2 (∂T ) , =c ≤c

T ∈T

(h2T ||R||2L2 (T )

)

hT ||J||2L2 (γ) )

(6)

T ∈T

where hT denote the diameter of T , apart from the constant c all quantities on the righthand side of (6) can be calculated explicitly. Then, an H 1 -type local error estimator ηT,H 1 associated with the element T ∈ T is defined by [6, 10] 1 1 ηT,H 1 = (h2T ||R||2L2 (T ) + hT ||J||2L2 (∂T ) ) 2 . 2

(7)

Also, duality arguments can be used to derive L2 -type a posteriori error estimators; here we just give the expression, details can be found in [10] and the references cited therein. The L2 -type local error estimator ηT,L2 is defined by 1

ηT,L2 = (h4T ||R||2L2 (T ) + h3T ||J||2L2 (∂T ) ) 2 ,

(8)

which is similar to (7); the only difference is a higher-order scaling in the mesh size; this reflects the expectation of a high-order rate of convergence with respect to the L2 norm. 3. Adaptive spatial discretization In an adaptive procedure, a good error estimator and an appropriate refinement strategy are two important issues needed to be considered. In this work, the H 1 -type and L2 -type local error estimators defined in Eq. (7) and (8) are used and shown to be asymptotically exact by numerical experiments. Then 7

the adaptive refinement strategy using the BLMG method at each refinement level will be discussed in this section, and firstly the mesh size function needs to be determined based on the error estimators. 3.1. Determination of the size function (l)

(l)

Let T (l) denote the triangulation of Ω with vertices {xi }N i=1 at the refine(l)

(l)

ment level l, N (l) is the number of vertices. Let ηT,H 1 and ηT,L2 represent the

corresponding local H 1 -type and L2 -type error estimators on T ∈ T (l) defined

in Eq. (7) and (8), respectively. According to the two types of error estimators, the size functions of the current triangular mesh will be modified and the new (l+1)

size functions for the l + 1 refinement level are denoted by hH 1 Let the patch mon vertex

(l) xi

(l) T˜i

(l+1)

and hL2

.

represent the union of the elements that share the com-

of T (l) , and the element patch is indicated by our local mesh

generation method which will be discussed in Section 3.3. Then the a poste(l)

riori error estimate associated with the vertex xi weighted average (l)

ηi =

P

is defined by area/volume

(l)

(l) T ∈T˜i

ηT |T |

P

(l) T ∈T˜i

|T |

(9)

,

(l+1)

where |T | denote the area/volume of element T . The new size functions hH 1 (l+1)

and hL2

are respectively determined by a piecewise linear function on Ω such (l)

that for any vertex xi of T (l) , (l)

(l+1)

(l)

hH 1 (xi ) =

(l)

hH 1 (xi ) (l)

f (ηi,H 1 )

(l)

(l+1)

and hL2

(l)

(xi ) =

(l)

hL2 (xi ) (l)

f (ηi,L2 )

(10)

,

(l)

(l)

where h(l) (xi ) is the average of the size of the elements sharing the vertex xi (l)

(l)

at the lth refinement level. ηi,H 1 and ηi,L2 denote the error estimates at the (l)

vertex xi when respectively using the H 1 -type and L2 -type a posteriori error (l)

(l)

estimators. For simplicity, both of them are denoted by ηi , and f (ηi ) is a user function defined by (l)

(l)

f (ηi ) = max( min( 8

ηi , a), b), c(l)¯ η (l)

(11)

where η¯(l) denote the mean of the a posteriori error estimates over all the elements at the lth level, c(l) is a coefficient generally close to 1.0, and the parameters a and b denote the mesh refinement and coarsening factor thresholds, respectively, i.e. a > 1, 0 < b < 1. Next, from Eq. (10) and (11), the value of f will be discussed from the following three kinds of cases: 1. If

(l)

ηi c(l)¯ η (l)

≥ a, then f = a holds, this implies that the mesh size in regions

with large errors is refined to be h(l+1) = a1 h(l) . 2. If b ≤

(l)

ηi c(l)¯ η (l)

< a , we can have f =

(l)

ηi , c(l)¯ η (l)

then h(l+1) =

c(l)¯ η (l) (l) h (l) ηi

can

be derived, and the new mesh size changes continuously according to the error.

3. If

(l)

ηi
to be h

, then f = b holds, this implies that the mesh size is enlarged

= 1b h(l) .

From the above three kinds of cases, the inequalities b ≤ f ≤ a can be

obtained, which indicate that the mesh size at one refinement level is respectively reduced or increased to be h(l+1) = a1 h(l) or h(l+1) = 1b h(l) at most. In order to study how to choose the values of parameters a and b, we consider the Poisson equation (1) on the unit square Ω = [0, 1]2 with the exact solution [33] u = 0.0005x2(1 − x)2 y 2 (1 − y)2 e10x

2

+10y

,

which exhibits an exponential peak near the corner point (1, 1). Since b is the mesh coarsening factor threshold, and it is usually chosen from 0.5 to 1.0. Suppose b=0.7 and let a=1.5, 2.0, 3.0 and 4.0, respectively, the error L2 norm versus the number of nodes in the adaptive process is demonstrated in the left of Fig. 1. From the left of Fig. 1, it can be seen that the meshes will deliver a small error when a → 2.0 in the case of the same number of nodes.

Similarly, when a=1.5, the smaller error norm can be obtained relative to the case of a=2.0. However, it will require more iterations to achieve almost the

same L2 error level when a=1.5 compared to the case of a = 2.0, which can be seen in Table 1. This demonstrates that it is reasonable for a to be 2.0 in practice. Then let a=2.0 and b is chosen to be 0.6, 0.7 and 0.8, respectively. The 9

−1

−2.5

a=2.0 b=0.7 a=3.0 b=0.7 a=4.0 b=0.7

−1.5

a=2.0 b=0.6 a=2.0 b=0.7 a=2.0 b=0.8

−3 log10(||.||L2)

log10(||.||L2)

−2

−2.5

−3

−3.5

−3.5

−4

1

1.5

2

2.5 log10N

3

3.5

−4

4

2

2.2

2.4

2.6

2.8 3 log10N

3.2

3.4

3.6

3.8

Figure 1: The influence of parameters a and b on the convergence rate of the finite element solution. Left: Fix b and change a. Right: Fix a and change b.

Table 1: The iterations when achieving almost the same L2 error level

Adaptive refinement using hH 1 N

l

||e||L2

a = 1.5, b = 0.7

2805

11

0.00012

a = 2.0, b = 0.7

2351

7

0.00016

a = 3.0, b = 0.7

3100

6

0.00014

a = 4.0, b = 0.7

4520

6

0.00013

error L2 norm versus the number of nodes in the adaptive process is shown in the right of Fig. 1. From the right of Fig. 1, we can see that the error L2 norm of the finite element solution is much smaller when b = 0.7 and b =0.8 relative to the case of b = 0.6. Thus it is relatively appropriate to set b = 0.7 or b =0.8. In the following numerical examples, the mesh coarsening factor threshold b is chosen to be 0.7. Meanwhile, the ratio

(l)

ηi η ¯(l)

between the error at the vertex and the average

error reflects the distribution of errors. In order to accelerate the convergence rate of the finite element solution, the status of the error distribution can be artificially enlarged through division of the ratio

(l)

ηi η ¯(l)

by the coefficient c(l),

where c(l) < 1.0, such that the triangles in the vicinity of the vertex with large 10

error are refined quickly. (l+1)

Here, we will refer to the mesh size functions hH 1

(l+1)

and hL2

as the H 1 -

based and L2 -based size functions, respectively. The parameters a and b are chosen to be 2.0 and 0.7 based on the above analysis. From Eq. (10) and (11), it can be observed that if f = 2.0, then h(l+1) = (l) refine the mesh near the vertex xi . If (l) the triangles around the vertex xi will

1 (l) 2h

holds, and we should

f = 0.7 holds, h(l+1) =

1 (l) 0.7 h ,

then

be coarsened. This coincide with the

principle of coarsening the meshes in regions where errors are relatively small and refining the meshes in regions where the errors are relatively large. In addition, since the size function for any point x ∈ Ω is defined by linear

interpolation with respect to the current mesh T (l) , when calculating the size

function value of the point x, it is required to search its nearest neighbor and find out the element where the point x lies, and the operation of searching is very time consuming generally. In fact, this work is efficiently accomplished with the software package “ANN” [34] which is based on kd-tree search strategies. If the point x is located on the vertex of T (l) , the size function value is obtained

based on Eq. (10). If the point x is located on an edge or in an element of T (l) ,

the size function value is calculated by linear interpolation. 3.2. Node placement method

For the mesh-based numerical methods, the accuracy and convergence properties are dependent on the size and shape of the elements, which consequently rely on the distribution of the corresponding nodes set. In the adaptive process, the node placement method by bubble simulation and its improved strategies are adopted to optimize an existing node distribution (refer to Refs. [26–28] for details), such that the refined/coarsened nodes set satisfies the requirement of the new size function. The flowchart is given and shown in Fig. 2, and the main procedures are discussed respectively in the following. (l)

The current mesh vertices {xi }N i=1 are considered as the centres of the pack-

ing bubbles, N (l) is the number of bubbles at the refinement level l, and the radius of the bubble centered at xi is

h(xi ) 2 .

11

Then according to Newton’s second

The current mesh vertices are considered as initial bubbles

Compute new positions of bubbles

Boundary constraints

N

Add/ delete bubbles

Renew adjacency list every 5 steps

Is the overlap ratio of bubble is satisfied?

N

Y

Whether the specified iterations reaches? Y

The centers of bubbles are used as nodes

Figure 2: The flowchart of the node placement process.

law of motion, the motion equation of each bubble is m¨ xi + c˜x˙ i = f˜i ,

i = 1, 2, · · · , N (l)

(12)

where m is the mass of bubble, xi is the center of bubble i, and f˜i is the resultant force exerted on bubble i by its neighbors. Additionally, −˜ cx˙ i is a viscous

damping force from the system, which makes the bubble system converge to

a stable configuration. Further, c˜ is the viscosity coefficient which gradually increases instead of being taken as a constant, which helps to speed up the convergence of the bubble system [28]. For any two neighboring bubbles, the interaction force can be approximated by [21]:   k 1.25w3 − 2.375w2 + 1.125
0 ≤ w ≤ 1.5 0 f˜ (w) =  0 1.5 < w

(13)

where w is the ratio of the real distance and the desired distance between two

bubbles, and the desired distance for two bubbles i and j is 12

h(xi ) 2

+

h(xj ) 2 .

To solve the equations (12), a second-order Euler predictor-corrector method is used here, and an initial velocity and an initial position are required for the bubbles. The initial position is the current position of bubbles, and the initial velocity is set to be zero. The zero initial velocity has an impact on the simulation process, but has no effect on the final distribution of bubbles. In addition, the boundary nodes are generated firstly and fixed during the next dynamic simulation, the interior nodes which conform with the boundary nodes in some sense are then determined. For the interior nodes which move outside the domain during dynamic simulation, they must be projected into the domain and kept a certain distance from the boundaries. Due to the short-range interaction force between bubbles, an adjacency list for the bubble i including its neighboring bubbles located within the cutoff radius r = 1.5σ is defined, where σ is the ideal distance between bubble i and its neighbours. When calculating the resultant force, only the interaction forces exerted by adjacent bubbles from the adjacency list need to be considered. And each adjacency list is usually updated every 5 steps. Moreover, the number of bubbles is adjusted by deleting the bubbles, the overlap ratios of which are too large, or adding new bubbles near the bubbles with overlap ratios that are too small during dynamic simulation. The overlap ratio [21] for each bubble is defined as αi =

1 Xn (2ri + rj − lij ) j=0 ri

where ri and rj are the ideal radii of bubbles i and j, respectively, ri = and rj =

h(xj ) 2 .

(14) h(xi ) 2

lij is the real distance between the centers of bubbles i and

j, and n is the size of the adjacency list associated with bubble i. In fact, the overlap ratio αi indicates the number of adjacent bubbles. In the ideal case, the standard overlap ratios of bubbles on a line, on a surface or in the internal volume are 2, 6 and 12, respectively. Through the setting of threshold values of the overlap ratios for the bubbles, the overlap ratio of bubbles are computed per a fixed iteration. Then the bubbles with large overlap ratios will be deleted, while the ones with small overlap ratios will have new bubbles added in their 13

Central node P

Satellite element Satellite nodes

Figure 3: A local mesh T˜P associated to a central node P

neighborhood. Finally, a node set meeting the requirements of the size function is obtained. 3.3. Bubble-type local mesh generation(BLMG) Some notations are firstly introduced in this section. In a Delaunay triangulation, the elements which share the common node P compose an element patch denoted by T˜P . The common node P is designated as a central node, and each element in the patch T˜P is called a satellite element of the central node P , which can be seen in Fig. 3. A node of a satellite element, other than the central node is referred to as a satellite node. The node placement procedure mentioned above provides not only a highquality node set, but also the adjacency list information which make the mesh generation very simple and efficient. The adjacency list related to each node has the following two features [29]: 1. The adjacency list of each node includes all of its satellite nodes in the corresponding Delaunay triangulation of the node set; 2. The adjacency list contains a small number of non-satellite nodes (no more than 2). The first feature implies that the satellite nodes of each central node can be found out in a small range. The second feature indicates that the satellite nodes of each central node can be obtained only through deleting a handful of non-satellite nodes from the corresponding adjacency list, and then the local mesh associated with the central node is generated. Next, taking the central 14

Central node P

Pj+1

Adjacent nodes Other nodes Pj

Pj-1

Lines to be detected Test line

Figure 4: Intersecting edges caused by connecting central node and its non-satellite nodes

node P as an example, a simple strategy of removing its non-satellite nodes is given as follows. 1. Connecting the central node P with all the adjacent nodes from its corresponding adjacency list; 2. Sorting the nodes from the adjacency list in counterclockwise order,we can get the sequence · · · Pj−1 , Pj , Pj+1 · · · which is shown in Fig. 4;

3. Checking the node Pj (j = 1, 2, · · · , n) whether a satellite node of the central node P , n is the size of the adjacency list.

If the nodes Pj−1 and Pj+1 don’t lie in each other’s adjacency lists, then the node Pj is a satellite node of the central node P . Otherwise, according to step 1, the nodes Pj−1 and Pj+1 are connected and the edge Pj−1 Pj+1 intersects with the edge P Pj , then Delaunay criteria can be used to check the position relationship between the node Pj and the circumscribed circle of ∆P Pj−1 Pj+1 . If Pj lies outside the circumscribed circle of ∆P Pj−1 Pj+1 , then Pj is a non-satellite node of the central node P , deleting Pj from the adjacency list of the node P , otherwise, Pj is a satellite node of P . After deleting the non-satellite nodes of each central node, only the satellite nodes are left in the adjacency list, and the local mesh associated with the central node P , i.e. the patch T˜P is obtained. In fact, the union of the patch T˜P is a Delaunay triangulation of the analysis domain. A simple example of local mesh generation is given in Fig. 5. The initial mesh through connecting all nodes with their adjacent nodes is shown in the left of Fig. 5. After deleting the 15

−0.5

−0.5

−1

−1

−1.5

−1.5

−2

−2

−2.5

−2

−1.5

−1

−0.5

−2.5

−2

−1.5

−1

−0.5

Figure 5: Local meshes generation example. Left: Connecting all nodes with their related adjacent nodes. Right: After deleting non-satellite nodes

non-satellite nodes of each central node, the resulting mesh is obtained which is shown in the right of Fig. 5. From Fig. 5 and the deletion strategy of non-satellite nodes, it can be seen that the BLMG method is very simple to implement, and a high-quality nodes set with a simple way of connecting nodes ensure the quality of triangle elements. 4. BLMG-based adaptive finite element methods Let Ω denote the given domain, Nl denote the number of nodes at the refinement level l, and Nmax > 0 is the allowable maximal number of mesh vertices. The flowchart of the adaptive finite element method based on BLMG is given as follows: 1. Generate an initial coarse triangulation T (0) of Ω and set l = 0; N0 is the number of vertices of T (0) .

2. Solve Eq. (1) by finite element method (FEM) on T (l) . If Nl > Nmax or l is more than the given refinement level, terminate; otherwise, go to step 3. (l)

(l)

3. Compute the local error estimators ηT,H 1 (or ηT,L2 ) for each element T ∈

T (l) . Determine the new size function of each mesh vertex using Eqs.

16

(l+1)

(9-11), then the meshsize field hH 1 interpolation with respect to T

(l)

(l+1)

(or hL2

) is obtained by linear

. (l+1)

4. Optimal node placement according to the new size function hH 1

(or

(l+1) hL2 ).

5. Build the satellite element patch corresponding to each central node with the BLMG method, and let T (l+1) represents the union of the patches. Set l = l + 1, and then go to step 2.

5. Numerical examples In this section, several numerical tests are considered to illustrate the effectiveness of the BLMG-based AFEM. Using the two-type error estimators mentioned in Section 2 in the adaptive process, the comparative study about the number of nodes, mesh quality, finite element solution accuracy and convergence rate during iteration are carried out. In Eq. (12) and (13), we set the coefficient k0 = 1.0, the mass coefficient m = 1.0 and the viscous coefficient c˜ = 3.8429, respectively. Initial coarse meshes are chosen to be produced using the “TRIANGLE” package [35] and are subsequently refined by the BLMG-based AFEM repeatedly. In addition, the exact solutions are pre-determined in all the experiments; they contain different forms of singularities that are frequently encountered in practical applications. The convergence rate CR with respect to the norm || · || at the refinement level l is roughly computed by [10, 12] CR =

2 log(||eh,l ||/||eh,l−1 ||) , log(Nl−1 /Nl )

where Nl denotes the number of nodes and eh,l denotes the error u − uh at the

refinement level l. Let ||eh ||L2 (Ω) and |eh |H 1 (Ω) denote the L2 norm and semiH 1 norm. Let ||er ||L2 denote the relative error in the L2 norm, i.e. ||er ||L2 = ||u−uh ||L2 ||u||L2

.

We apply the commonly used q measure [10, 36] to evaluate the geometric

quality of triangular meshes where for any triangle T , q is defined to be twice 17

−3

x 10 3 2

Error

1 0 −1 −2 −3 1 0.8

1 0.6

0.8 0.6

0.4

0.4

0.2 Y

0.2 0

0

X

−3

x 10 3 2

Error

1 0 −1 −2 −3 1 1 0.8

0.5

0.6 0.4

Y

0

0.2 0

X

−3

x 10 3 2

Error

1 0 −1 −2 −3 1 1 0.8

0.5

0.6 0.4

Y

0

0.2 0

X

Figure 6: Adaptive BLMG meshes for Example 1. Left: 431, 869 and 3453 nodes using the size function hH 1 . Right: The distribution of absolute error at the corresponding refinement levels.

the ratio of the radius RT of the largest inscribed circle and the radius rT of the smallest circumscribed circle, i.e. q(T ) =

(a1 + a2 − a3 )(a2 + a3 − a1 )(a1 + a3 − a2 ) 2RT = , rT a1 a2 a3

where a1 , a2 and a3 are side lengths of T . 5.1. Example 1: With a shock wave The first example is to solve Eq. (1) on domain Ω = [0, 1]2 . The right-hand side f and the Dirichlet boundary conditions are chosen such that the exact 18

solution is [37] u=

1 1+

e20(x+y)−25

.

The solution function represents a shock wave along the line y = 1.25 − x. The

BLMG-based AFEM is used to obtain a convergent solution.

The initial coarse mesh is chosen to be uniform mesh with h(0) = 0.1 (l=0), (l+1)

then the adaptive meshes generated by the BLMG method using hH 1 (l+1) hL 2

and

are shown in the left of Fig. 6 and Fig. 7, respectively, where the three

meshes corresponding to the 2nd, 3rd and 5th levels of adaptive refinement are shown, and the absolute error distributions corresponding to the 2nd, 3rd and 5th adaptive meshes are also shown in the right of Fig. 6 and Fig. 7, respectively, (l+1)

where it can be clearly seen that the absolute error for using both hH 1 (l+1) hL 2

and

is strictly reduced and almost equally distributed on the elements at last.

From the left of Fig. 6 and Fig. 7, we can also observe that the BLMG meshes (l+1)

generated using the size function hH 1

tend to distribute the nodes in a slightly (l+1)

less uniform manner than those generated using the spacing function hL2 can be explained by the fact that

(l+1) hH 1

has larger variations than does

, this

(l+1) hL 2 .

Table 2 gives the information about mesh quality, solution errors and convergence rates at all refinement levels with different refinement strategies for Example 1. From Table 2, it can be seen that the BLMG-based adaptive meshes (l+1)

remain well shaped at all refinement levels for both size functions hH 1 (l+1) hL 2 ,

and

which is supported by the values of qmin and qavg . From the 3rd refine-

ment level, the average qualities of the meshes are greater than 0.96; this is a desirable property when solving PDEs with the finite element method. Meanwhile, the corresponding plots of the error norm (||eh ||L2 (Ω) and |eh |H 1 (Ω) )

versus the number of nodes are shown in Fig. 8. As expected, the adaptive refinement strategies lead to a smaller error than that obtained with the uniform refinement method. From Table 2 and Fig. 8, the adaptive method using (l+1)

hL 2

as the size function generates approximate solutions uh having smaller (l+1)

||eh ||L2 (Ω) relative to those obtained using hH 1 , this also has been verified in

the following three examples. The convergence rates asymptotically tend to be

19

−3

x 10 6 4

Error

2 0 −2 −4 −6 1 1 0.8

0.5

0.6 0.4

Y

0

0.2 0

X

−3

x 10 1.5 1

Error

0.5 0 −0.5 −1 −1.5 −2 1 1 0.8

0.5

0.6 0.4 0.2

0

Y

0

X

−3

x 10 1.5 1

Error

0.5 0 −0.5 −1 −1.5 −2 1 1 0.8

0.5

0.6 0.4

Y

0

0.2 0

X

Figure 7: Adaptive BLMG meshes for Example 1. Left: 428, 832 and 3298 nodes using the size function hL2 . Right: Absolute error distribution at the corresponding refinement levels.

2 and 1 for ||eh ||L2 (Ω) and |eh |H 1 (Ω) , respectively, which implies the reliability

and efficiency of our adaptive procedure.

In Table 3, one can observe that for the same PDE, and the same exact solution, in [37], the L2 norm of error with the three error estimators are respectively 1.4 × 10−3 , 3.5 × 10−4 and 3.4 × 10−4 on the adaptive anisotropic

mesh with 684, 693 and 714 vertices. However, from Table 1, for the BLMG-

based adaptive method, we have that ||eh ||L2 (Ω) = 3.738 × 10−4 on the refined mesh with 428 vertices; this means that we can achieve higher accuracy with less degree of freedom compared to the method in [37].

20

Table 2: Mesh quality, solution errors, and convergence rates for different refinement strategies for Example 1.

l

N

qmin

||e||L2

qavg

|e|H 1

CR

CR

||er ||L2

Adaptive refinement using hH 1 0

126

0.710

0.931

6.712e-03

0.0

1.129e-01

0.0

8.205e-03

1

205

0.538

0.938

7.780e-04

8.855

3.833e-02

4.437

9.472e-04

2

431

0.404

0.944

4.374e-04

1.550

2.258e-02

1.425

5.322e-04

3

869

0.494

0.961

1.979e-04

2.262

1.468e-02

1.229

2.406e-04

4

1728

0.435

0.969

9.853e-05

2.029

9.399e-03

1.296

1.296e-04

5

3453

0.450

0.975

5.700e-05

1.581

6.824e-03

0.925

8.710e-05

6

6873

0.484

0.982

2.593e-05

2.289

4.228e-03

1.391

4.671e-05

Adaptive refinement using hL2 0

126

0.710

0.931

6.712e-03

0.0

1.129e-01

0.0

8.205e-03

1

202

0.679

0.944

7.929e-04

9.051

3.744e-02

4.675

9.655e-04

2

428

0.448

0.945

3.738e-04

2.003

2.852e-02

0.725

4.547e-04

3

832

0.494

0.966

1.436e-04

2.879

1.547e-02

1.840

1.745e-04

4

1724

0.650

0.976

7.315e-05

1.851

8.821e-03

1.543

8.722e-05

5

3298

0.508

0.980

3.466e-05

2.303

5.928e-03

1.225

4.230e-05

6

6821

0.488

0.981

1.917e-05

1.630

4.508e-03

0.754

2.006e-05

5.2. Example 2: With two point singularities Consider the problem   −∇ · (∇u) + 10.0u = f in Ω ,  u = g on ∂Ω

the domain Ω = [−1, 1]2 . The exact solution u assumes the following form [10]: u(x, y) =

1.0 2

2

(x − 0.5) + (y − 0.5) + 0.01

−

1.0 2

2

(x + 0.5) + (y + 0.5) + 0.01

,

the exact solution u given above is a smooth function which reaches its maximal values at two points (0.5,0.5) and (-0.5,-0.5), but decreases rapidly away from 21

−2

−0.8 Uniform refinement Adaptive refinement using h

2

Adaptive refinement using h

1

−1

L

−2.5

Uniform refinement Adaptive refinement using h

2

Adaptive refinement using h

1

L

H

H

−1.2 −1.4 log10(|.|H1)

log10(||.||L2)

−3

−3.5

−4

−1.6 −1.8 −2 −2.2

−4.5 −2.4 −5

2

2.5

3

3.5

4

−2.6

4.5

2

2.5

3

log N

3.5

4

4.5

log N

10

10

Figure 8: Error norms versus the number of nodes for different refinement strategies for Example 1. Left: ||e||L2 . Right: |e|H 1 . Table 3: Comparison of the error L2 norm between the method in [37] and the BLMG-based AFEM.

Adaptive results in [37] Edge-based er-

Least

square

ror estimator

Hessian recov-

Full error esti-

The present method

mator

ery N ||e||L2

684

693

714

428

1.4 × 10−3

3.5 × 10−4

3.4 × 10−4

3.7 × 10−4

these maximums and thus has large gradients close to the two points. Applying the BLMG-based AFEM, the total 8 adaptive meshes are generated for the convergence of the whole procedure. Three meshes at different (l+1)

refinement levels are chosen for both size functions hH 1

(l+1)

and hL2

respec-

tively. See Fig. 9, which clearly shows the successive refinement around the two singular points. Also, the error distributions of elements for different refinement (l+1)

levels are respectively given in Fig. 10 (using the size function hH 1 ) and Fig. (l+1)

11 (using the size function hL2

), where it can be clearly seen that the element

errors are asymptotically equal-distributed on the elements. As in Example 1, the adaptive meshes optimized with the BLMG method

22

Figure 9: Adaptive BLMG meshes for Example 2. Top: 561, 5040 and 10319 nodes using the size function hH 1 . Bottom: 632, 5188 and 10559 nodes using the size function hL2 .

1

0.2

0.2

0.5 0

0

−0.5

Error

Error

Error

0

−0.2

−0.2

−1 −1.5 1

−0.4 1 0.5

−0.4 1 0.5

1 0.5

0 −0.5 −1

−1

1 0.5

0

0

−0.5

−0.5 −1

Y

X

0.5

0

−0.5

−0.5 −1

Y

1 0.5

0

0

−0.5 −1

Y

X

−1

X

Figure 10: Plots of error distribution on the mesh with 561, 2460 and 5040 nodes using hH 1 .

1

0

−0.5

0.4 0.2 0 Error

0.2

0 Error

Error

0.4 0.5

−0.2 −0.4

−1

−0.4

−0.6

−1.5 1

−0.6

−0.8 1 0.5

1 0.5

0

0

−0.5 Y

−0.8 1 0.5

1

−1

X

0.5

0

0

−0.5

−0.5 −1

−0.2

Y

0.5

−1

1

X

0.5

0

0

−0.5

−0.5 −1

Y

−0.5 −1

−1

X

Figure 11: Plots of error distribution on the mesh with 632, 1298 and 2654 nodes using hL2 .

23

Table 4: Mesh quality, solution errors, and convergence rates for different refinement strategies for Example 2.

l

N

qmin

qavg

||e||L2

CR

|e|H 1

CR

||er ||L2

Adaptive refinement using hH 1 0

147

0.688

0.954

9.785

0.0

60.411

0.0

0.475

1

196

0.557

0.924

7.053e-01

18.283

17.080

8.782

3.013e-02

2

313

0.211

0.890

3.672e-01

2.790

6.815

3.926

1.523e-02

3

561

0.550

0.952

2.883e-01

0.829

2.971

2.845

8.017e-03

4

1178

0.539

0.966

1.659e-01

1.489

1.602

1.666

5.231e-03

5

2460

0.662

0.977

7.647e-02

2.104

0.985

1.320

2.080e-03

6

5040

0.499

0.983

3.640e-02

2.070

0.623

1.278

1.205e-03

7

10319

0.475

0.985

1.903e-02

2.179

0.425

1.069

6.122e-04

Adaptive refinement using hL2 0

147

0.688

0.954

9.785

0.0

60.411

0.0

0.475

1

204

0.527

0.912

8.760e-01

14.729

9.614

11.218

3.744e-02

2

310

0.447

0.922

3.217e-01

4.788

7.221

1.368

1.336e-02

3

632

0.6570

0.966

1.690e-01

1.808

2.886

2.575

6.967e-03

4

1298

0.451

0.973

7.826e-02

2.139

1.975

1.054

3.228e-03

5

2654

0.490

0.978

4.343e-02

1.647

1.194

1.407

1.790e-03

6

5188

0.500

0.982

1.958e-02

2.377

0.799

1.197

8.071e-04

7

10559

0.455

0.984

1.103e-02

1.614

0.538

1.112

4.548e-04

24

1

2 Uniform refinement Adaptive refinement using h

2

Adaptive refinement using h

1

L

0.5

Uniform refinement Adaptive refinement using h

2

Adaptive refinement using h

1

L

H

H

1.5

log10(|.|H1)

log10(||.||L2)

0

−0.5

1

0.5

−1 0

−1.5

−2

2

2.5

3

3.5

4

−0.5

4.5

2

2.5

log N

3

3.5

4

4.5

log N

10

10

Figure 12: Error norms versus the number of nodes for different refinement strategies for Example 2. Left: ||e||L2 . Right: |e|H 1 .

have nice qualities, which are demonstrated by the value of qmin and qavg contained in Table 4. Table 4 also includes information about solution errors, and convergence rates at all refinement levels for different refinement strategies for Example 2. The quality of the mesh is always very good for both size func(l+1)

tion hH 1

(l+1)

and hL2

, even though the mesh sizes vary greatly in the domain;

the corresponding plots of the error norms (||eh ||L2 (Ω) and |eh |H 1 (Ω) ) versus the

number of nodes are given in Fig. 12. From Fig. 12, we can observe that the BLMG-based adaptive methods are much more efficient relative to the uni(l+1)

form refinement strategy, and the BLMG-based adaptive method using hL2

as the size function generates approximate solutions having smaller ||eh ||L2 (Ω) (l+1)

relative to those obtained using hH 1 ; on the other hand, the latter generates approximate solutions with smaller |eh |H 1 (Ω) . 5.3. Example 3: With layer singularity In this example, in order to test the capability of the BLMG-based AFEM to handle boundary layer singularities, the same PDE as in Example 1 is solved. The exact solution u is chosen to be [12] u(x, y) =

x2

1.0 1.0 + 2 . + 0.01 y + 0.01

25

Figure 13: Adaptive BLMG meshes for Example 3. Top: 1008, 4534 and 8810 nodes using the size function hH 1 . Bottom: 1072, 4840 and 9417 nodes using the size function hL2 .

Let Ω = Ω1 ∩ Ω2 − Ω3 , where Ω1 = [0, 4]2 , Ω2 = {(x, y)|x2 + y 2 ≤ 17} and Ω3 =

{(x, y)|(x − 2)2 +(y − 2)2 < 0.49} so that Ω is a nonconvex region which includes

two boundary layers along the x-axis and the y-axis. Fig. 13 displays the refined

meshes at some levels generated by the BLMG-based adaptive method using the (l+1)

size functions hH 1

(l+1)

and hL2

, respectively. All triangles are well shaped at all

refinement levels, an observation that is supported by the values of qmin and qavg listed in Table 5. Plots of the error norm (||eh ||L2 (Ω) and |eh |H 1 (Ω) ) versus the number of

nodes for different refinement strategies are presented in Fig. 14. From Table

5 and Fig. 14, similar to that in the first two experiments, the same conclusion about the convergence rate can be obtained. Also, the BLMG-based adap(l+1)

tive method using the size function hL2

generates approximate solutions with (l+1)

smaller ||eh ||L2 (Ω) relative to those obtained using hH 1 . On the other hand, the former and the latter generates the approximate solutions with almost the same values |eh |H 1 (Ω) .

26

Table 5: Mesh quality, solution errors, and convergence rates for different refinement strategies for Example 3.

l

N

qmin

qavg

||e||L2

CR

|e|H 1

CR

||er ||L2

Adaptive refinement using hH 1 0

384

0.583

0.932

20.567

0.0

88.580

0.0

0.244

1

487

0.529

0.929

1.430

22.437

15.641

14.595

1.774e-02

2

1008

0.440

0.943

0.354

3.837

13.519

0.401

4.326e-03

3

2133

0.485

0.956

0.192

1.637

7.463

1.585

2.330e-03

4

4534

0.474

0.957

0.100

1.720

4.652

1.254

1.219e-03

5

8810

0.467

0.971

0.046

2.329

2.936

1.386

5.624e-04

6

17613

0.504

0.968

0.021

2.290

1.876

1.293

2.775e-04

Adaptive refinement using hL2 0

384

0.583

0.932

20.567

0.0

88.580

0.0

0.244

1

511

0.527

0.926

1.4828

18.408

18.014

11.149

1.839e-02

2

1072

0.512

0.953

0.292

4.391

11.257

1.269

3.561e-03

3

2332

0.362

0.959

0.116

2.379

7.075

1.195

1.407e-03

4

4840

0.497

0.9628

0.0518

2.205

4.802

1.062

6.292e-04

5

9417

0.441

0.9763

0.025

2.137

2.948

1.466

3.115e-04

6

18802

0.526

0.980

0.013

2.029

1.954

1.189

1.603e-04

27

1.5

2 Uniform refinement Adaptive refinement using h

2

Adaptive refinement using h

1

1.8

L

1

Uniform refinement Adaptive refinement using h

2

Adaptive refinement using h

1

L

H

H

1.6 0.5

log (|.| 1)

H

1.2

10

log10(||.||L2)

1.4 0

−0.5

1 0.8

−1 0.6 −1.5

−2

0.4

2.6

2.8

3

3.2

3.4 log N

3.6

3.8

4

4.2

0.2

4.4

2.6

2.8

10

3

3.2

3.4 log N

3.6

3.8

4

4.2

4.4

10

Figure 14: Error norms versus the number of nodes for different refinement strategies for Example 3. Left: ||e||L2 . Right: |e|H 1 .

5.4. Example 4: L-shape domain problem In this section we compare the performance of the BLMG-based adaptive strategy with a highly optimized h-adaptive software bamg [38], where the bisection technique and mesh smoothing are used. For this purpose, we consider the Laplace equation −∆u = 0 in the L-shape domain Ω = [−1, 1] ×

[−1, 1]\([−1, 0) × [−1, 0)). The homogeneous Dirichlet boundary condition is imposed on the whole boundary ∂Ω so that the exact solution is [18] π 2 u(r, θ) = r2/3 sin (θ + ), 3 2 where (r, θ) is a polar representation of a point in the L-shape domain. Derivatives of the solution are singular at origin (0, 0) ( see, e.g., Ref. [39]). For illustration, several steps of the BLMG-based adaptive algorithm are shown in Fig. 15. It can be easily seen that the meshes near the singularity at the re-entrant corner get significantly refined. A comparison of convergence rate of the relative error in L2 norm is shown in Fig. 16, where the x-axis represents the number of degrees of freedom (DOFs). Here, since we use the linear element approximation for all numerical examples in this study, the number of degrees of freedom is also the node population. From Fig. 16, we note that the asymptotic convergence rate for the adaptive methods is much better than that in the uniform refinement case, both the bamg and BLMG-based adaptive refinements 28

Figure 15: Refined adaptive meshes generated by the BLMG-based adaptive method. Top, from left to right: T (0) with 23 nodes and T (3) with 200 nodes. Bottom, from left to right: T (6) with 2023 nodes and the enlarged view near the singular point.

get good approximate solutions as h → 0. However, to obtain similar accuracy, i.e. the relative error ||er ||L2 = 2.0e − 05, it uses 7987 DOFs in the case of bamg

software based on bisection technique (see Level 10 in Table 6), while it only uses 4135 DOFs in the adaptive BLMG case (see Level 7 in Table 6). Thus, the BLMG-based adaptive refinements can obtain much better approximate solution than the bamg at almost the same number of degrees of freedom. 5.5. Time consumption In the BLMG-based adaptive method, well-shaped triangle meshes at all refinement levels can be obtained, and the mesh refining and coarsening can also be acquired simultaneously in the same framework, which can be seen in the left of Figs. 6 and 7. In addition, from Tables 2 and 4-6, we can observe that the number of nodes is increased proportionally with the refinement level. Accordingly, the CPU time spent on generating adaptive meshes is increased linearly with respect to the refinement level, which can be seen in Fig. 17.

29

Table 6: Comparison of the relative errors in L2 norm using two adaptive strategies for Example 4.

bamg using bisection technique

BLMG-based adaptive method

l

DOF s

||er ||L2

l

DOF s

||er ||L2

0

21

1.964e-02

0

23

1.436e-02

2

47

6.506e-03

1

49

7.398e-03

4

137

1.787e-03

2

97

2.861e-03

5

258

8.678e-04

3

200

1.149e-03

6

507

4.169e-04

4

426

3.333e-04

8

2006

8.941e-05

5

923

1.280e-04

9

4001

4.573e-05

6

2023

5.004e-05

10

7987

2.118e-05

7

4135

2.001e-05

−1

10

Uniform refinement Mesh adapting software bamg Adaptive BLMG method −2

2

Relative error in L norm

10

−3

10

−4

10

−5

10

1

10

2

10

3

10 Number of DOF

4

10

5

10

Figure 16: Convergence history for the L-shape domain problem.

30

700

450 400

600

350

500 CPU time(s)

CPU time(s)

300 250 200

400

300

150

200 100

100 50 0

0 0

200

400 600 800 1000 The number of nodes per level

1200

1400

0

500

1000 1500 The number of nodes per level

2000

2500

Figure 17: CPU time spent on the BLMG-based adaptive method. Left: The point singularity problem in Fig. 1. Right: L-shape domain problem.

In general, the time consumption for the presented method is mainly spent on the four aspects: mesh size modification, the node placement procedure, the node-based local mesh generation and the solution of linear equations. Typically in the last example of the L-shape domain, the time for equilibrium of bubbles is several times greater than the solution time. Actually it largely depends on the details of coding, and could be further reduced if the coding is implemented in an efficient manner. 6. Conclusions In this paper, we present an adaptive mesh refining algorithm for elliptic PDEs that combines a posteriori error estimation with the bubble-type local mesh generation. The BLMG-based AFEM assures that errors are very well equidistributed over the triangles, and at all levels of refinement, the triangles remain very well shaped. Several numerical examples are demonstrated to verify that the proposed BLMG-based adaptive finite element method is robust and effective for the numerical solutions of elliptic equations, and optimal convergence rates with respect to both H 1 and L2 norms are obtained. Though the initial study on the application of BLMG-based AFEM to the numerical solutions of PDEs carried out here is for a classical model of two31

dimensional linear second order elliptic equation, we expect greater advantages when BLMG is applied to solve more complex systems of equations (such as high order equations). Also, since the satellite elements are generated independently in the local region of each node, this node-based local mesh generation method has potential parallelism, and the study of the parallel BLMG-based adaptive finite element method will be our future work. Acknowledgements This research was supported by National Natural Science Foundation of China (No: 11471262), and the Doctorate Foundation of Northwestern Polytechnical University. We would like to thank the referees for their many valuable comments and suggestions, which have led to a significantly improved presentation of this paper. References [1] R. Stevenson, An optimal adaptive finite element method, SIAM Journal on Numerical Analysis 42 (5) (2005) 2188–2217. [2] P. Binev, W. Dahmen, R. Devore, Adaptive finite element methods with convergence rate, Numerische Mathematik 97 (2) (2004) 219–268. [3] X. Xu, H. Chen, R. H. Hoppe, Optimality of local multilevel methods on adaptively refined meshes for elliptic boundary value problems, Journal of Numerical Mathematics 18 (1) (2010) 59–90. [4] I. Babuvˇska, W. Rheinboldt, Error estimates for adaptive finite element computations, SIAM Journal on Numerical Analysis 15 (4) (1978) 736–754. [5] R. Bank, A. Weiser, Some a posteriori error estimators for elliptic partial differential equations, Mathematics of Computation 44 (170) (1985) 283– 301.

32

[6] T. Gr¨atsch, K. J. Bathe, A posteriori error estimation techniques in practical finite element analysis, Computers and Structures 83 (4-5) (2005) 235–265. [7] K. Segeth, A review of some a posteriori error estimates for adaptive finite element methods, Mathematics and Computers in Simulation 80 (8) (2010) 1589–1600. [8] M. Ainsworth, J. Oden, A posteriori error estimation in finite element analysis, John Wiley & Sons, New York, 2000. [9] R. Bank, R. Smith, Mesh smoothing using a posteriori error estimates, SIAM Journal on Numerical Analysis 34 (3) (1997) 979–997. [10] L. L. Ju, M. Gunzburger, W. D. Zhao, Adaptive finite element methods for elliptic PDEs based on conforming centroidal voronoi-Delaunay triangulations, SIAM Journal on Science Computing 28 (6) (2006) 2023–2053. [11] G. Nicolas, T. Fouquet, Adaptive mesh refinement for conformal hexahedral meshes, Finite Elements in Analysis and Design 67 (2013) 1–12. [12] Y. Q. Huang, H. F. Qin, D. S. Wang, Q. Du, Convergent adaptive finite element method based on centroidal voronoi tessellations and superconvergence, Communications in Computational Physics 10 (2) (2011) 339–370. [13] T. Nguyen-Thoi, G.R.Liu, H. Nguyen-Xuan, C. Nguyen-Tran, Adaptive analysis using the node-based smoothed finite element method (NS-FEM), International Journal for Numerical Methods in Biomedical Engineering 27 (2) (2011) 198–218. [14] R. Stevenson, The completion of locally refined simplicial partitions created by bisection, Mathematics of Computation 77 (261) (2008) 227–241. [15] M. Paszynski, D. Pardo, V. M. Calo, A direct solver with reutilization of lu factorizations for h-adaptive finite element grids with point singularities, Computers & Mathematics with Applications 65 (8) (2013) 1140–1151. 33

[16] J. Helsing, R. Ojala, Corner singularities for elliptic problems: Integral equations, graded meshes, quadrature, and compressed inverse preconditioning, Journal of Computational Physics 227 (20) (2008) 8820 – 8840. [17] L. Demkowicz, W. Rachowicz, P. Devloo, A fully automatic hp-adaptivity, Journal of Scientific Computing 17 (1-4) (2002) 117–142. [18] L. Demkowicz, Computing with hp-adaptive finite elements: Volume 1 one and two dimensional elliptic and maxwell problems, CRC Press, 2006. [19] J. Maciej, Z. Grzegorz, On some hp-adaptive finite element method for natural vibrations, Computers and Mathematics with Applications 66 (11) (2013) 2376–2399. [20] J. Lang, W. M. Cao, W. Z. Huang, R. D. Russell, A two-dimensional moving finite element method with local refinement based on a posteriori error estimates, Applied Numerical Mathematics 46 (1) (2003) 75–94. [21] K. Shimada, D. Gossard, Automatic triangular mesh generation of trimmed parametric surfaces for finite element analysis, Computer Aided Geometric Design 15 (3) (1998) 199–222. [22] O. C. Zienkiewicz, J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, International Journal for Numerical Methods in Engineering 24 (2) (1987) 337–357. [23] J. Kim, H. Kim, B. Lee, S. Im,

Adaptive mesh generation by bubble

packing method, Structural Engineering and Mechanics 15 (1) (2003) 135– 149. [24] S. Chung, S. Kim, A remeshing algorithm based on bubble packing method and its application to large deformation problems, Finite Elements in Analysis and Design 39 (4) (2003) 301–324. [25] O. C. Zienkiewicz, J. Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. part 1: The recovery technique, International Journal for Numerical Methods in Engineering 33 (7) (1992) 1331–1364. 34

[26] Y. Liu, Y. F. Nie, W. W. Zhang, L. Wang, Node placement method by bubble simulation and its application, CMES-Computer Modeling in Engineering & Sciences 55 (1) (2010) 89–109. [27] Y. F. Nie, W. W. Zhang, N. Qi, Y. Li, Parallel node placement method by bubble simulation, Computer Physics Communications 185 (3) (2014) 798–808. [28] N. Qi, Y. F. Nie, W. W. Zhang, Acceleration strategies based on an improved bubble packing method, Communications in Computational Physics 16 (1) (2014) 115–135. [29] W. W. Chen, Y. F. Nie, W. W. Zhang, L. Wang, A fast local mesh generation method about high-quality node set, Chinese Journal of Computational Mechanics 29 (5) (2012) 704–709. [30] W. Zhang, Y. Nie, Y. Li, A new anisotropic local meshing method and its application in parametric surface triangulation, Computer Modeling in Engineering & Sciences(CMES) 88 (6) (2012) 507–529. [31] W. Zhang, Y. Nie, L. Wang, Bubble meshing method for two-parametric surface, Chinese Journal of Computational Physics 29 (1) (2012) 43–50. [32] J. Ku, A. H. Schatz, Local a posteriori estimates on a nonconvex polygonal domain, SIAM Journal on Numerical Analysis 50 (2) (2012) 906–924. [33] K. Y. Kim, Guaranteed a posteriori error estimator for mixed finite element methods of elliptic problems, Applied Mathematics and Computation 218 (24) (2012) 11820–11831. [34] S. Arya, D. M. Mount, N. S. Netanyahu, R. Silverman, A. Y. Wu, An optimal algorithm for approximate nearest neighbor searching fixed dimensions, Journal of the ACM 45 (6) (1998) 891–923. [35] J. Shewchuk, Triangle: a two-Dimensional quality mesh generator and Delaunay triangulator http://www.cs.cmu.edu/quake/triangle.html. 35

[36] P. Persson, G. Strang, A simple mesh generator in MATLAB, SIAM REVIEW 46 (2) (2004) 329–345. [37] W. Z. Huang, L. Kamenski, J. Lang, A new anisotropic mesh adaptation method based upon hierarchical a posteriori error estimates, Journal of Computational Physics 229 (6) (2010) 2179–2198. [38] F. Hecht, The mesh adapting software: bamg, INRIA report, 1998. [39] A. Szymczak, A. Paszy´ nska, P. Gurgul, M. Paszy´ nski, Graph grammar based direct solver for hp-adaptive finite element method with point singularities, Procedia Computer Science 18 (2013) 1594–1603.

36