finite element method

Chin-Tien Wu,

Zhilin Li

Ming-Chih Lai,

October 30, 2009

Abstract

In this paper, an adaptive mesh refinement technique is developted and analyzed for the non-conforming immersed finite element (IFE) method pro-posed in [25]. The IFE method was developed for solving the second order elliptic boundary value problem with interfaces across which the coefficient may be discontinuous. The IFE method was based on a triangulation that does not need to fit the interface. One of the key ideas of IFE method is to modify the basis functions so that the natural jump conditions are satis-fied across the interface. The IFE method has shown to be order of O(h²) and O(h) in L² norm and H¹ norm, respectively. In order to develop the adaptive mesh refinement technique, additional priori and posterior error es-timations are derived in this paper. Our new a priori error estimation shows that the generic constant is only linearly proportional to ratio of the diffusive coefficients β⁻ and β⁺, which improves the corresponding result in [25].

∗Corresponding author. Department of Applied Mathematics, National Chiao-Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan. [email protected]

†Center for Research in Scientific Computation & Department of Mathematics, North Carolina State University, Raleigh, NC. 27695-8205.

‡Center of mathematical modeling and Scientific computing, National Chiao-Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan.

1

We also show that a posteriori error estimate similar to the one obtained by Bernardi and Verf¨urth [4] holds for the IFE solutions. Numerical examples support our theoretical results and show that the adaptive mesh refinement strategy is effective for the IFE approximation.

1 Introduction

The main purpose of this paper is to develop adaptive mesh refinement techniques for the immersed finite element (IFE) method proposed in [25]. Along this line, we also discuss the a priori and a posteriori error estimation for the immersed finite element method. The IFE method was developed for the following interface problem:

−∇ · (β∇u) = f, (x, y) ∈ Ω

u |∂Ω = g, (1)

together with the natural jump conditions on the interface ˜Γ:

[u] |Γ˜= 0, (2)

[βu_n] |Γ˜= 0. (3)

Here, see the sketch in Fig.1, Ω ⊂ R²is a convex polygonal domain, the interface Γ is a curve separating Ω into two sub-domains Ω˜ ⁻, Ω⁺such that Ω = Ω⁻∪Ω⁺∪˜Γ, and the coefficient β(x, y) is a piecewise constant function defined by

β(x, y) =

½ β⁻, (x, y) ∈ Ω⁻, β⁺, (x, y) ∈ Ω⁺.

The interface problem considered here appears in many engineering and sci-ence applications. The immersed finite element (IFE) space was first introduced in [25], in which some preliminary analysis and numerical results are reported, and has been shown its capability on handling interface problems with nonho-mogeneous interface jump conditions [with a nonzero constant value on the right hand of (2) and/or (3)] by either simply modifying the IFE space or reducing the interface problem to a new problem with homogeneous interface jump conditions

2

−3 −2 −1 0 1 2 3

−3

−2

−1 0 1 2 3

\tex{$\tilde{Gamma}$}

\tex{$\beta^+$}

\tex{$\beta^−$}

\tex{$\Omega$}

\tex{$\Omega^+$}

\tex{$\Omega^−$}

Γ˜ β⁺

β⁻

Ω Ω⁺

Ω⁻

Figure 1: A rectangular domain Ω = Ω⁺∪ Ω⁻with an immersed interface ˜Γ. The coefficients β(x) may have a jump across the interface.

via the usual homogenization technique based on a change of variable [23]. Some related work can be found in [?, 18, 19, 26].

The basic idea of the immersed finite elements is to form a partition =_h in-dependent of interface ˜Γ so that partitions with simple regular structures can be used to solve an interface problem with a rather complicated or varying interface.

Obviously, triangles in a partition can be separated into two classes:

• Non-interface triangles: The interface ˜Γ either does not intersect with this triangle, or it intersects with this triangle but does not separate its interior into two nontrivial subsets.

• Interface triangles: The interface ˜Γ cuts through its interior.

In a non-interface triangle, the standard linear polynomials is employed as local basis functions. However, in an interface triangle, a piecewise linear polynomial is defined in the two subsets formed by the interface in a way that the functions satisfy the natural jump conditions (either exactly or approximately) on the in-terface and retain specified values at the vertices of the inin-terface triangle. The immersed finite element space defined over the whole domain Ω can then be con-structed through the standard finite element assembling procedure. We refer the readers to [?, 9–11, 14, 17, 22, 24] for more background materials about immersed interface and immersed finite element methods as well as their applications.

3

Without loss of generality, we assume that the triangles in the partition have the following features:

(H1): If ˜Γ meets one edge of a triangle at more than two points, then the edge is part of ˜Γ.

(H2): If ˜Γ meets a triangle at two points, then these two points must be on different edges for this triangle.

In order to obtain error estimates, we assume that the underlying mesh is fine enough such that the interface can be approximated by a line segment with a small perturbation in a magnitude of O(h²). Furthermore, the source function f and the interface ˜Γ are assumed to be smooth enough such that the weak solution of the problem (1) can be approximated by a piecewise C²function. These requirements lead to our third hypothesis:

(H3): The segment of the interface ˜Γ in a triangle T ∈ =h is defined by a piecewise C²function and the function space C²(T ) is dense in H²(T ).

It is well known that the standard finite element method (FE) with linear finite elements can be used to solve such elliptic interface problems [see [3, 5, 6] and the references therein]. However, in order to achieve the optimal O(h²) accuracy in the numerical solutions, an interface fitted grid is needed. In applications with nontrivial interfaces or the time-varying interfaces, this restriction prevents the Galerkin method with linear finite elements from working efficiently since mesh moving or re-meshing is required. On the other hand, although the mesh moving and re-meshing may produce extra technical difficulties and computation over-head for the standard FE method, the standard FE method has a great advantage on increasing the accuracy of the numerical solutions at low cost through the adaptive mesh refinement process. In the adaptive mesh refinement process, first an error indicator ηT used to pin point the locations with large error is computed on each element in a given triangulation. Second, the elements in which the error indicator has large value are marked for refinement according to a given marking strategy. A heuristic marking strategy is the maximum marking strategy where an element T^∗ will be marked for refinement if η_T^∗ > θ max_{T ∈=h}η_T, with a prescribed threshold 0 ≤ θ ≤ 1. Some other marking strategies can also be seen in [13]. Finally, the marked triangles are divided into sub-triangles by rules such as the regular

4

refinement algorithm or the longest side bisection algorithm [15] [16]. An ap-proximate solution is then computed on the refined mesh. The above procedure can be repeatedly applied until the accuracy of the approximated solution is sat-isfied. The theoretical foundation of the mesh refinement strategy is based on the a posteriori error estimation proposed by Babuˇska and Rheinboldt [1] and further developed by many researchers such as Zienkiewicz [27], Bank and Weiser [2], and Verf¨urth [20, 21]. The convergence of the adaptive mesh refinement process has been shown by Morin, Nochetto and Siebert [12].

It has been shown that the IFE interpolation errors on a uniform fixed (such as Cartesian) partition is of the order of O(h) in the H¹ norm and of the order of O(h²) in the L^∞ and L² norms under the hypothesis (H₁), (H₂) and (H₃) [26].

In this work, we obtain the same order of the error estimations and further show that the generic constants in these estimations are linearly proportional to the ratio max

n ρ,¹_ρ

o

of the diffusion coefficients, here ρ = ^β_β⁻+. The a posteriori esti-mations of the finite element solutions mentioned above are obtained mostly on fitted grids. Recently, A. Hansbo and P. Hansbo propose an unfitted finite element method for the elliptic interface problem. The same order of a priori error esti-mations is obtained and an a posteriori estimator is proposed [8]. Here, we also derive an a posteriori error estimation for the IFE method based on the methodol-ogy developed by Verf¨urth [4]. Our numerical results support the effectiveness of the proposed a posteriori error estimation.

This paper is organized as follows. In section 2, we show the existence and uniqueness of the element IFE basis function and derive some auxiliary inequal-ities that are needed for the error estimation in section 3. We derive the a priori error estimations and the a posteriori error estimation in section 3 and present our numerical results in section 4. Finally, we draw our conclusions in section 5.