A partially penalty immersed Crouzeix-Raviart finite element method for interface problems.

The elliptic equations with discontinuous coefficients are often used to describe the problems of the multiple materials or fluids with different densities or conductivities or diffusivities. In this paper we develop a partially penalty immersed finite element (PIFE) method on triangular grids for anisotropic flow models, in which the diffusion coefficient is a piecewise definite-positive matrix. The standard linear Crouzeix-Raviart type finite element space is used on non-interface elements and the piecewise linear Crouzeix-Raviart type immersed finite element (IFE) space is constructed on interface elements. The piecewise linear functions satisfying the interface jump conditions are uniquely determined by the integral averages on the edges as degrees of freedom. The PIFE scheme is given based on the symmetric, nonsymmetric or incomplete interior penalty discontinuous Galerkin formulation. The solvability of the method is proved and the optimal error estimates in the energy norm are obtained. Numerical experiments are presented to confirm our theoretical analysis and show that the newly developed PIFE method has optimal-order convergence in the [Formula: see text] norm as well. In addition, numerical examples also indicate that this method is valid for both the isotropic and the anisotropic elliptic interface problems.

Introduction

The elliptic equations with discontinuous coefficients are often used to describe phenomena appearing in material sciences and fluid dynamics when there are two or more distinct materials or fluids with different densities or conductivities or diffusivities. Since the solutions of these interface problems are required to satisfy interface jump conditions from conservation laws, it is difficult to find the exact solutions and construct high accuracy numerical methods. In addition, if the interface is smooth enough, then the solution of the interface problem is also smooth in individual regions where the coefficient is smooth. But due to the jump of the coefficient along the interface, the global regularity of the solution is usually low and the solution belongs to , . Therefore, it is difficult to achieve high accuracy by using standard finite element methods. Several articles are devoted to developing methods to solve these problems, such as fitted finite element methods [1-3], proposed by Chen and Zou in [4], discontinuous Galerkin methods [5, 6], weak Galerkin finite element methods [7, 8], proposed by Wang and Ye in [9-11], immersed boundary methods [12, 13], the extended/generalized finite element methods [14], immersed interface methods, which are a kind of finite difference methods, proposed by LeVeque and Li in [15], and immersed finite volume methods [16, 17].

In recent years, the immersed finite element (IFE) methods have been studied and found to be very effective for solving elliptic interface problems. The method was proposed on the uniform Cartesian triangular grids in [18, 19], according to the above mentioned immersed interface methods. Approximation capabilities of the nonconforming IFE spaces were studied in [20] and the convergence analysis of the IFE solutions was developed in [21, 22]. But we realize that it is very difficult to achieve an optimal-order -norm error estimate due to the strong nonconformity. So, for anisotropic elliptic interface problems, we presented the partially penalty IFE method in [23] by adding two penalty terms on the common edges of the adjacent interface elements to restrict the function jumps, and then deriving the optimal-order error estimates. The numerical experiments verify our theoretical results. Moreover, this idea of partially penalization can also be seen in [24], which is on rectangular grids for solving isotropic elliptic interface problems. And there are other developed forms [25-27]. In addition, the penalty idea exists in other unfitted methods, too. In [28], a ghost penalty is added in unfitted finite element methods to recover the condition number of the stiffness matrix. And an unfitted method based on the symmetric interior penalty discontinuous Galerkin method was proposed to discretize elliptic interface problems in [29]. Moreover, there are other forms of IFE methods, such as the symmetric and consistent IFE method [30] and the augmented IFE method [31, 32].

Recently, Kwak et al. [33] proposed an immersed finite element method based on piecewise linear Crouzeix-Raviart type polynomials on a uniform triangular grid. In this reference, the authors use the integral averages on the edges as degrees of freedom to weaken their nonconformity on the common edge of two adjacent interface elements. Zhang, in his thesis [34], proposed the nonconforming rotated IFE method based on rectangular grids by using both integral averages and midpoint values on the edges as degrees of freedom.

However, the optimal convergence is still not very easily obtained. And inspired by the work of Ji et al. [28], a penalty on interface edges is added to ensure that the condition number of the resulting linear system is independent of how the interface cuts through the mesh. In this paper, we develop the partially penalty immersed finite element (PIFE) method with the Crouzeix-Raviart type polynomial spaces to solve the anisotropic flow models in which the diffusion coefficient is a piecewise definite-positive matrix, although we consider isotropic elliptic problems with piecewise scalar coefficients in major length for simplicity. For the anisotropic elliptic interface problems, the construction of the Crouzeix-Raviart type IFE spaces is given in our paper [35]. We prove that the piecewise Crouzeix-Raviart type polynomials satisfying the jump conditions on two types of interface elements can be uniquely determined by integral averages on the edges as degrees of freedom. Then we construct the Crouzeix-Raviart type IFE spaces on interface elements and give the partially penalty immersed Crouzeix-Raviart finite element schemes. We prove the solvability of the method and obtain its optimal convergence analysis. Numerical experiments show that our method is valid not only for isotropic elliptic interface problems but also for anisotropic elliptic interface problems.

The rest of the paper is organized as follows. In the next section, some preliminaries and notations are introduced. In Section 3, we construct IFE spaces based on the Crouzeix-Raviart elements on triangular grids and investigate properties of these nonconforming IFE spaces. In Section 4, we are devoted to defining the PIFE method based on the Crouzeix-Raviart type IFE spaces and prove the solvability of the method. Then we carry out the optimal convergence analysis. In Section 6, some numerical experiments are performed to indicate the optimal-order convergence of our PIFE method for both isotropic and anisotropic elliptic interface problems. Finally, we conclude in the last section.

Preliminaries and notations

Consider the second order elliptic interface problem: where is a discontinuous coefficient, , . The convex polygonal domain consists of and , ; see Figure 1 (cited from [23]) for illustration. Assume that the interface is smooth enough () and the solution u satisfies the following jump conditions across the interface Γ: where the jump value and denotes the unit normal vector of Γ pointing from to .

Figure 1

The geometry of the computational domain Ω.

For the sake of simplicity and convenience, in the following analysis, we assume is a piecewise constant function defined by where , .

Of course, the coefficient can also be a symmetric definite-positive matrix as and the corresponding analysis can be performed similarly. In the last section, we also conduct numerical examples for this case.

Let be the usual regular triangulation of the domain Ω. and denote the collection of interface elements and the collection of the non-interface elements, respectively. It is also called a non-interface element if Γ intersects with this triangle but does not separate its interior into two nontrivial subsets. Assume that the interface meets the edges of an interface element at no more than two intersections. Such an assumption is reasonable if the step size is sufficiently small.

Let be the collection of all edges in the triangulation . and denote the interior edges and boundary edges, respectively. Moreover, the sets of the interface edges and non-interface edges are denoted by and . Obviously, here we have and also .

In the following, we define the jump and average values of a function u on the edges. For every interior edge , where and are the two elements sharing the common edge e and the unit normal vector of e is assumed to point from to . For every boundary edge , where is the element such that is one of its edges. Usually we omit the subscript in and if there is no confusion.

For the analysis, we introduce the following spaces on the whole domain Ω: and for any , where and are the usual Sobolev spaces.

For every interface element , we also introduce the following spaces: and for any , where is the norm of , .

The Crouzeix-Raviart type IFE space and its properties

In this section, we introduce the local Crouzeix-Raviart type basis functions for both the non-interface elements and the interface elements, define the IFE spaces over the whole domain, and then investigate properties of the nonconforming Crouzeix-Raviart type IFE spaces.

Construction of immersed Crouzeix-Raviart finite element spaces

In this subsection, we develop the nonconforming IFE spaces with integral averages on the edges as degrees of freedom. To make sure that the flux jump conditions can be weakly enforced on the smooth interface, we derive the IFE functions on usual elements instead of the reference elements used in [33]. The standard nonconforming linear Crouzeix-Raviart functions are used on non-interface elements. On interface elements, these functions are locally modified to satisfy the jump conditions.

For a non-interface element , we simply use the standard linear Crouzeix-Raviart type polynomials as local basis functions, and use to denote spaces spanned by the three basis functions on T, where , are three edges of T and δ is the Kronecker function.

For the space on non-interface element , we have the well-known approximation property as follows [36, 37]: where is the interpolation operator defined by

To construct the local immersed Crouzeix-Raviart finite element spaces on interface elements, we consider a typical triangle , in which the three vertices are , and . We assume that the interface curve Γ intersects T at two different points D and E. And the segment separates T into two subsets and with .

There are two types of interface elements depending on the location of the intersection points D and E. We call an element T a Type-1 interface element if Γ intersects with T at two square edges, or a Type-2 interface element if Γ intersects with T at a square edge and a bevel edge. See Figure 2 for an illustration of the two different types of interface elements, where, without loss of generality, we assume that intersection points D, E satisfy for a Type-1 interface element with , and for a Type-2 interface element with , .

Figure 2

Type-1 (left) and Type-2 (right) interface triangular elements.

There is a small region in these two types of interface elements, whose area is of order since the interface is a curve and the interface is perturbed in a magnitude of [4]. Therefore, such a perturbation will only affect the solution and the interpolation function to an order of , which will not impact on the convergence accuracy of the method whose approximation spaces are selected as piecewise linear polynomials in this paper.

On each of these interface triangular elements, , for given values , , the piecewise linear function ϕ can be defined by satisfying where , denote integral averages of ϕ on the edges , i.e., , and is a unit normal vector to .

As for the piecewise linear function ϕ defined above, we have the following results:

Theorem 3.1

For arbitrary diffusion coefficient β and interface location, the piecewise linear function on the interface element defined by (3.3) and (3.4) is uniquely determined by the given three values , , .

Proof

We shall prove the theorem on Type-1 and Type-2 interface elements. In both cases, it is proved by involving , , , as unknowns into a linear system and then showing its coefficient matrix is non-singular for arbitrary interface location, reflected by a and b.

Firstly, we consider the Type-1 interface elements, where , with . By the definition of and (3.4a) we obtain, through a simple calculation, the first three equations with respect to unknowns,

The continuity conditions (3.4b) and (3.4c) imply the fourth and fifth equations,

The last equation is derived from the flux continuity condition (3.4d) along , where we used and .

We summarize these six equations and reformulate them into the following matrix form: where , and the coefficient matrix

In the following we will show that the coefficient matrix is non-singular or its determinant is non-zero. Following a tedious calculation, we have By , we have . Then implies that for this case.

For the Type-2 interface elements, where , with , , by the definition of and (3.4a), we derive the first three equations similarly to the derivation of (3.5),

The continuity conditions (3.4b) and (3.4c) imply two equations as below, The last equation is obtained from the flux continuity condition (3.4d), where we used and .

We summarize these six equations and write in the matrix form: where , and the coefficient matrix

In the following we will prove that the coefficient matrix is non-singular. Following a tedious calculation, we have

To simplify the proof, we let with , then is rewritten as

By , we have . Also by , we derive and . Therefore implies in this case, which completes the proof. □

Remark 3.2

Theorem 3.1 tells us that on each interface element, the piecewise linear function defined by (3.3) and (3.4) is continuous across and uniquely determined by its integral averages on the edges as degrees of freedom.

Then, for an interface element , we can define the local Crouzeix-Raviart type IFE space by where denotes the inverse matrix of that is derived in the proof of Theorem 3.1.

Remark 3.3

Theorem 3.1 also holds if the coefficient is a symmetric definite-positive matrix defined by (2.3) and its proof is provided in another paper [35].

Finally, we conclude this subsection by defining the global immersed Crouzeix-Raviart type finite element spaces and on the whole domain Ω as

Properties of the Crouzeix-Raviart type IFE spaces

In this subsection, we present several properties for the local IFE space on the interface element and the global IFE space on the domain Ω.

Although for Crouzeix-Raviart type IFE functions in the flux jump condition is enforced on the line segment , they actually satisfy a weak flux jump condition along the actual interface curve , which is stated in the following lemma.

Lemma 3.4

For an interface triangle , each linear function satisfies the flux jump condition on in the following weak sense:

Let ϕ be a function in . By Green’s formula we have Combining the equality above with the flux continuity of ϕ along we reach the conclusion. □

It is well known that the trace inequalities are important for the finite element analysis. So here we present the trace inequalities stated on p.23 of [38], Noting that , the above two inequalities both hold on the non-interface element . However, is a subspace of only, not a subspace of . Therefore, we need to establish a trace inequality similar to (3.21) for the Crouzeix-Raviart type IFE functions on interface elements.

Theorem 3.5

There exists a constant C depending on the diffusion coefficient β only, such that on any , where e is an interface edge of T, and n is the unit outer normal vector of T.

To prove Theorem 3.5, which follows the idea of proof for Theorem 3.3 in [34], we need the lemma below.

Lemma 3.6

There exist two constants , , which are both dependent on the discontinuous coefficient β but independent of the mesh size h and the interface location (reflected by a, b), such that, for every function ϕ defined by (3.3), we have

We shall prove the lemma on Type-1 and Type-2 interface elements, respectively.

Firstly, we consider the Type-1 interface elements, where , with . According to equations (3.6) and (3.7), we have By a tedious calculation, we derive Then, by using and , we get the inequality for , And in the same way, the inequality for is derived: That is to say, we have

If we express , , by , , by use of (3.6) and (3.7), then we get the inequalities for similarly:

Therefore, we derive the conclusion (3.23) with

Next we consider the Type-2 interface elements, where , with , .

According to equations (3.11) and (3.12), we have

By a tedious calculation, we obtain where .

Using and , we derive the inequality for , Then the inequality for is derived similarly: That is to say, we have

If we express , , by , , by applying (3.11) and (3.12), then we can get the inequalities for similarly, Therefore, we also derive the conclusion for Type-2 interface elements. □

For any interface element , we can divide it into four congruent triangles , by connecting the three midpoints of three edges of T. There must exist a subset which is completely inside of either or , and we denote it by T̃.

If , noting that , we have Applying the left of the estimate (3.23) in the above lemma, we derive So we have the estimate

If , by the right of the estimate (3.23), we obtain

Therefore, no matter which part of T the subset T̃ belongs to, we have

Now we give the proof for Theorem 3.5.

Without loss of generality, we consider with and . Let , by the above result (3.24), we obtain which completes the proof. □

For any , we let be such that where , are three edges of T. We call the interpolant of u in , and we have its error estimate which is stated in [33], where is the norm in defined by (2.5).

We can naturally extend the interpolation operator such that and .

Then, by (3.1) and the above estimation (3.25), we obtain the following result.

Theorem 3.7

For any , there exists a constant such that where .

The Crouzeix-Raviart PIFE method

We multiply the elliptic equation (2.1a) by a test function , integrate over each element , and apply the Green’s formula, Here, is the unit outward normal of T.

Summing the above equation over all elements, we obtain Let denote the unit normal vector of e, then Applying the algebraic identity leads to Then substituting the second term in (4.1) by the above two equalities, we obtain

Assume that u is smooth enough, implies that and thus the normal flux is continuous across every interior edge e. Hence, . Therefore the third term in the above equation is zero and the above equation becomes where we have used the definition of the jump and average values for functions on boundary edges and the relationship .

In the following, similar to [23], we give the partially penalty immersed finite element method with Crouzeix-Raviart elements and prove the existence and uniqueness of its discrete solution. In the derivation of the PIFE formulation, we have followed the idea of the interior penalty discontinuous Galerkin scheme [38]. The meaning of partially penalty is that we add penalty terms only to the interface edges in order to restrict the jump of IFE functions across their interface edges.

Since u is continuous in the interior of Ω, on each interior edge e. We add two penalty terms whose value is zero for the exact solution u defined only on interface edges to (4.2), and obtain

We assume that on non-interface edges, the quantity is not very large which suggests to ignore this term in our scheme. Then we have where we have used the relationship .

Now we define the PIFE formulation: find such that where In the above formulation, the penalty parameter , the choice of the power will be discussed later, and the parameter ε may take the value −1, 0, or 1.

Here, we also call the scheme (4.4) by the same name as reference [23] (shown in Remark 4.1), although the degrees of freedom used are different.

Remark 4.1

If , the scheme (4.4) is called the symmetric PIFE method; if , the scheme (4.4) is called the nonsymmetric PIFE method; if , the scheme (4.4) is called the incomplete PIFE method.

Subsequently, we shall prove the existence and uniqueness of the solution to (4.4). Define the energy norm on by It is easy to check that it is indeed a norm on .

In the following theorem, the coercivity of the bilinear form with respect to the energy norm is stated, which is similar to that in [23] although the degrees of freedom used are different.

Theorem 4.2

The coercivity of the bilinear form :

is coercive on , that is, , .

and are coercive if and if is bounded below by a constant that depends on the diffusion coefficient β and the constant in the trace inequality.

Proof

Note that the coercivity result is trivial for and the corresponding coercivity constant is . Indeed, Hence, we focus on the other two cases or 0 below.

For every interface edge shared by the interface elements and , using Cauchy-Schwarz’s inequality twice and the trace inequality stated in Theorem 3.5, we derive where and are the maximum edge length of and , respectively, and the common edge e satisfies , .

If satisfies the condition and we assume, without loss of generality, that . Then, from the above inequality, we obtain Summing the above (4.6) over all the interface edges, and using Cauchy-Schwarz’s inequality once more and the Young inequality, for , we get Thus, we obtain a lower bound for , where C is the constant in the trace inequality stated in Theorem 3.5.

Choosing δ and in the above inequality such that then we have the coercivity result for the cases or 0 as follows: For instance, if and if and choosing large enough (for example, if and if ), then we have the coercivity result with .

Summarizing the results above, we complete the proof. □

Now we are ready to derive the existence and uniqueness of the PIFE solution as we have done in the paper [23] (which takes values at three vertices as degrees of freedom).

Theorem 4.3

Assume that (I) or (II) holds true: Then the PIFE solution of the discrete problem (4.4) exists and is unique.

in the nonsymmetric () PIFE case, for all e;

in the symmetric () or incomplete () PIFE case, and is bounded below by a large constant for all e.

Since (4.4) is a linear problem in finite dimension, existence is equivalent to uniqueness. We only have to prove the uniqueness of the PIFE solution.

Assuming that there are two solutions and , the difference satisfies

Taking , by the coercivity results stated in Theorem 4.2, we have

This implies since is a norm on , which completes the proof. □

The convergence analysis of the PIFE method

In this section, we conduct the error estimate for the PIFE solution in the energy norm . And we firstly give the interpolant error estimate in the energy norm.

Lemma 5.1

For , there exists a constant C independent of h and interface location, such that

where we used that is bounded above for any when it is fixed.

By the interpolant error estimates given in (3.1) and (3.25) as well as the trace inequality (3.20), we obtain which completes the proof. □

Remark 5.2

The parameter in the energy norm (4.5) should satisfy in order to guarantee that the Crouzeix-Raviart type IFE space has the optimal approximation property in the energy norm.

Next, we come to derive the error estimate for the PIFE solution in the energy norm. In the beginning, we give the following lemma which is cited from Lemma 3 in [37], and we shall need it in the proof of the energy norm estimate.

Lemma 5.3

Let e be an edge of T. Then there exists a constant such that for all where .

Theorem 5.4

Let and be the solutions of the interface problem (2.1)-(2.2) and the discrete formulation (4.4), respectively. Assume , then there exists a constant independent of h and the location of the interface, such that the following optimal energy norm estimate holds:

The equation (4.3) implies that the solution u satisfies Subtracting (4.4) from the above equation, we obtain Let , then we have Choosing the test function and using the coercivity of the bilinear form in Theorem 4.2, we get Using Cauchy-Schwarz’s inequality and Young’s inequality, we have the estimate for , The approximation results (3.1) and (3.25) give the estimate for ,

Now we bound the term . By the Cauchy-Schwarz inequality, we have where and share the common edge e.

Using the trace inequality stated in Theorem 3.5 for the IFE function ξ, we have where .

And by the trace inequality (3.20) for the function η, we obtain Then combining with Young’s inequality and approximation result (3.25), we have

The term can simply be bounded using Cauchy-Schwarz’s and Young’s inequalities: Assume . Using the estimate (5.5) and the approximation result (3.25), we derive

Next, we come to bound .

From the definition of the Crouzeix-Raviart type IFE space in (3.18), for , we have which also implies for any edge e shared by and .

Note that , applying Lemma 5.3 by choosing , , , we obtain Therefore, applying the above estimate and (5.6), combining with the observation that is continuous across any interior edge , we obtain

Then by using Cauchy-Schwarz’s and Young’s inequalities, we derive

In order to conclude, it remains to bound the term . Similarly to what we did to the above terms, applying Cauchy-Schwarz’s inequality and Young’s inequality, we have Assume , then we obtain the estimate for , where we used , , which can be obtained similarly to Lemma 4.4 in [34] by using the approximation result (3.25).

Combining these five bounds above, we derive which implies the energy norm error estimate for ξ, that is, .

In addition, Remark 5.2 tells us that, for ,

Finally, we complete the proof by the triangle inequality  □

Remark 5.5

The error estimate in -seminorm can be derived as follows: that is to say,

To guarantee the optimal-order energy norm estimate (5.2) of the PIFE solutions, the parameter in our PIFE schemes (4.4) and the energy norm (4.5) should be chosen as . Hence, we assume for the numerical experiments in the next section.

Computational results

In this section, two kinds of numerical experiments are presented to illustrate the validity of our PIFE schemes with Crouzeix-Raviart elements. One kind is for the isotropic elliptic interface problems with piecewise constant coefficients, which we apply to verify our theoretical findings about the error estimate in the energy norm as well as to show the optimal convergence in norm. And the other is for the anisotropic elliptic interface problems in which the diffusion coefficients are piecewise symmetric definite-positive matrices. We conduct this kind of numerical experiments to show that the PIFE method proposed in this paper can also be applied to the problems with discontinuous tensor-coefficients, although we only analyze the problems with discontinuous scalar-coefficients for simplicity.

We consider the elliptic interface problem defined by (2.1a) and (2.2) with , Γ being a circle centered at origin with radius . Then the interface curve Γ separates Ω into two sub-domains and with

In our numerical experiments, we use the triangular grids, which are formed by partitioning Ω by congruent squares and then cutting every square along its two diagonals. Thus we derive our triangulation (illustrated in Figure 3) with mesh size , which is the maximum length of edges. The number of triangles in this partition is .

Figure 3

The triangular partition with .

Numerical results for the isotropic elliptic interface problem

In this subsection, we conduct the numerical experiment with piecewise constant-coefficients to verify our theoretical findings about the PIFE method in symmetric, nonsymmetric and incomplete forms defined in Remark 4.1.

Taking the example from [20], the boundary condition function and the source term are chosen such that the following function u is the exact solution:

In our computation, we select the penalty parameter in the discrete formulation as for the nonsymmetric PIFE scheme, for the incomplete PIFE scheme and for the symmetric PIFE scheme.

Then we present the numerical results for and solved by these three different PIFE schemes, which are reflected by different ε, in Table 1 and Table 2, respectively.

Table 1

Numerical results for with different PIFE schemes

Scheme	Grid	∥u−u_h∥_0,Ω	Order	|u−u_h|_1,h	Order
Nonsymmetric PIFE (ε = 1)	8 × 8	4.9337E − 3	-	8.0908E − 2	-
	16 × 16	1.1162E − 3	2.1441	4.3830E − 2	0.8844
	32 × 32	2.9009E − 4	1.9440	2.2917E − 2	0.9355
	64 × 64	6.5676E − 5	2.1431	1.1502E − 2	0.9946
	128 × 128	1.8919E − 5	1.7955	5.7744E − 3	0.9941
	256 × 256	4.4585E − 6	2.0852	2.8965E − 3	0.9954
Incomplete PIFE (ε = 0)	8 × 8	5.5911E − 3	-	8.1187E − 2	-
	16 × 16	1.3044E − 3	2.0997	4.4172E − 2	0.6087
	32 × 32	3.2578E − 4	2.0014	2.3182E − 2	0.9301
	64 × 64	6.9110E − 5	2.2369	1.1552E − 2	1.0048
	128 × 128	1.8967E − 5	1.8654	5.7761E − 3	1.0000
	256 × 256	4.4796E − 6	2.0820	2.8974E − 3	0.9953
Symmetric PIFE (ε = −1)	8 × 8	6.0212E − 3	-	8.1462E − 2	-
	16 × 16	1.3636E − 3	2.1426	4.4097E − 2	0.8854
	32 × 32	3.3744E − 4	2.0147	2.3284E − 2	0.9214
	64 × 64	7.1179E − 5	2.2451	1.1612E − 2	1.0037
	128 × 128	1.8993E − 5	1.9060	5.7769E − 3	1.0073
	256 × 256	4.4914E − 6	2.0802	2.8981E − 3	0.9952

Table 2

Numerical results for with different PIFE schemes

Scheme	Grid	∥u−u_h∥_0,Ω	Order	|u−u_h|_1,h	Order
Nonsymmetric PIFE (ε = 1)	8 × 8	1.0678E − 2	-	4.7925E − 1	-
	16 × 16	2.4087E − 3	2.1483	2.3967E − 1	0.9997
	32 × 32	5.7775E − 4	2.0597	1.1989E − 1	0.9993
	64 × 64	1.4172E − 4	2.0274	5.9947E − 2	1.0000
	128 × 128	3.5894E − 5	1.9812	2.9981E − 2	0.9996
	256 × 256	8.8877E − 6	2.0139	1.4993E − 2	0.9998
Incomplete PIFE (ε = 0)	8 × 8	1.0510E − 2	-	4.7886E − 1	-
	16 × 16	2.3468E − 3	2.1630	2.3937E − 1	1.0004
	32 × 32	5.7093E − 4	2.0393	1.1983E − 1	0.9982
	64 × 64	1.4131E − 4	2.0144	5.9947E − 2	0.9992
	128 × 128	3.5885E − 5	1.9775	2.9982E − 2	0.9996
	256 × 256	8.8873E − 6	2.0136	1.4993E − 2	0.9998
Symmetric PIFE (ε = −1)	8 × 8	1.0501E − 2	-	4.7883E − 1	-
	16 × 16	2.3443E − 3	2.1633	2.3935E − 1	1.0004
	32 × 32	5.7071E − 4	2.0383	1.1984E − 1	0.9981
	64 × 64	1.4128E − 4	2.0142	5.9948E − 2	0.9993
	128 × 128	3.5883E − 5	1.9772	2.9983E − 2	0.9996
	256 × 256	8.8910E − 6	2.0129	1.4993E − 2	0.9998

The data in Table 1 and Table 2 illustrate that, for each scheme, the PIFE solution converges to the exact solution with convergence order in the norm and in the -seminorm, which supports our theoretical findings stated in (5.7) and verifies the robustness and validity of this method.

Numerical results for the anisotropic elliptic interface problem

To show that the PIFE method proposed in this paper can also be applied to anisotropic elliptic interface problems, we provide a numerical example in this subsection, whose diffusion coefficients are piecewise symmetric definite-positive matrices , satisfying , , .

The boundary condition function and the source term are chosen such that the exact solution is which can easily be checked to satisfy the jump conditions (2.2).

In our computation, we select the parameter in the discrete formulation as for the nonsymmetric PIFE scheme, for the incomplete PIFE scheme and for the symmetric PIFE scheme.

Firstly, we present the numerical results for , and , solved by all the three PIFE schemes in Table 3 and Table 4, respectively. We can see the optimal-order accuracy in both the norm and the -seminorm, which implies that our PIFE method is valid for the anisotropic interface problem. In addition, we present the figures of the numerical solutions in these two numerical examples (the left one is for and the right one is for ), solved by the nonsymmetric () PIFE scheme under the triangular grid with (see Figure 4).

Table 3

Numerical results for ,

Scheme	Grid	∥u−u_h∥_0,Ω	Order	|u−u_h|_1,h	Order
Nonsymmetric PIFE (ε = 1)	8 × 8	1.2388E − 1	-	2.3454	-
	16 × 16	3.1375E − 2	1.9812	1.1815	0.9892
	32 × 32	7.9027E − 3	1.9892	5.9497E − 1	0.9898
	64 × 64	1.9817E − 3	1.9956	2.9851E − 1	0.9950
	128 × 128	4.9509E − 4	2.0010	1.4946E − 1	0.9980
	256 × 256	1.2382E − 4	1.9994	7.4772E − 2	0.9992
Incomplete PIFE (ε = 0)	8 × 8	1.2234E − 1	-	2.3267	-
	16 × 16	3.1249E − 2	1.9690	1.1786	0.9812
	32 × 32	7.8952E − 3	1.9848	5.9442E − 1	0.9875
	64 × 64	1.9834E − 3	1.9930	2.9841E − 1	0.9942
	128 × 128	4.9591E − 4	1.9998	1.4941E − 1	0.9980
	256 × 256	1.2407E − 4	1.9989	7.4762E − 2	0.9989
Symmetric PIFE (ε = −1)	8 × 8	1.2175E − 1	-	2.3184	-
	16 × 16	3.1205E − 2	1.9640	1.1778	0.9771
	32 × 32	7.8917E − 3	1.9834	5.9433E − 1	0.9867
	64 × 64	1.9836E − 3	1.9923	2.9840E − 1	0.9940
	128 × 128	4.9598E − 4	1.9997	1.4939E − 1	0.9981
	256 × 256	1.2411E − 4	1.9986	7.4760E − 2	0.9988

Table 4

Numerical results for ,

Scheme	Grid	∥u−u_h∥_0,Ω	Order	|u−u_h|_1,h	Order
Nonsymmetric PIFE (ε = 1)	8 × 8	1.8010E − 1	-	3.5464	-
	16 × 16	5.1864E − 2	1.7960	2.1268	0.7377
	32 × 32	1.4329E − 2	1.8558	1.1824	0.8469
	64 × 64	3.6979E − 3	1.9541	6.1435E − 1	0.9446
	128 × 128	9.5226E − 4	1.9573	3.1524E − 1	0.9626
	256 × 256	2.3939E − 4	1.9920	1.5897E − 1	0.9877
Incomplete PIFE (ε = 0)	8 × 8	1.7212E − 1	-	3.3852	-
	16 × 16	5.0836E − 2	1.7595	2.0585	0.7177
	32 × 32	1.4025E − 2	1.8578	1.1467	0.8440
	64 × 64	3.6820E − 3	1.9295	6.0852E − 1	0.9141
	128 × 128	9.5061E − 4	1.9536	3.1442E − 1	0.9526
	256 × 256	2.3906E − 4	1.9915	1.5868E − 1	0.9865
Symmetric PIFE (ε = −1)	8 × 8	1.7172E − 1	-	3.3561	-
	16 × 16	5.0887E − 2	1.7547	2.0499	0.7112
	32 × 32	1.4029E − 2	1.8589	1.1439	0.8416
	64 × 64	3.6839E − 3	1.9291	6.0784E − 1	0.9122
	128 × 128	9.5104E − 4	1.9537	3.1429E − 1	0.9516
	256 × 256	2.3912E − 4	1.9918	1.5865E − 1	0.9862

Figure 4

Numerical solutions under the triangular grid with .

Then we come to verify that our PIFE method can derive the same optimal-order accuracy even if the entries of the tensors () are negative. We conduct the numerical experiment for , and , and give the error results in Table 5 and Table 6, respectively.

Table 5

Numerical results for ,

Scheme	Grid	∥u−u_h∥_0,Ω	Order	|u−u_h|_1,h	Order
Nonsymmetric PIFE (ε = 1)	8 × 8	1.2304E − 1	-	2.3400	-
	16 × 16	3.1217E − 2	1.9788	1.1799	0.9878
	32 × 32	7.8667E − 3	1.9885	5.9443E − 1	0.9891
	64 × 64	1.9740E − 3	1.9946	2.9836E − 1	0.9944
	128 × 128	4.9317E − 4	2.0010	1.4940E − 1	0.9979
	256 × 256	1.2335E − 4	1.9993	7.4752E − 2	0.9990
Incomplete PIFE (ε = 0)	8 × 8	1.2054E − 1	-	2.3118	-
	16 × 16	3.1041E − 2	1.9572	1.1770	0.9739
	32 × 32	7.8555E − 3	1.9824	5.9417E − 1	0.9862
	64 × 64	1.9738E − 3	1.9927	2.9834E − 1	0.9939
	128 × 128	4.9332E − 4	2.0004	1.4937E − 1	0.9981
	256 × 256	1.2343E − 4	1.9988	7.4746E − 2	0.9988
Symmetric PIFE (ε = −1)	8 × 8	1.2010E − 1	-	2.3060	-
	16 × 16	3.1019E − 2	1.9530	1.1767	0.9707
	32 × 32	7.8549E − 3	1.9815	5.9414E − 1	0.9858
	64 × 64	1.9737E − 3	1.9927	2.9833E − 1	0.9939
	128 × 128	4.9330E − 4	2.0004	1.4936E − 1	0.9981
	256 × 256	1.2343E − 4	1.9988	7.4746E − 2	0.9987

Table 6

Numerical results for ,

Scheme	Grid	∥u−u_h∥_0,Ω	Order	|u−u_h|_1,h	Order
Nonsymmetric PIFE (ε = 1)	8 × 8	1.0816	-	23.4214	-
	16 × 16	3.9033E − 1	1.4704	16.6163	0.4952
	32 × 32	1.1628E − 1	1.7470	9.8040	0.7612
	64 × 64	3.1365E − 2	1.8904	5.3068	0.8855
	128 × 128	8.2005E − 3	1.9354	2.7690	0.9385
	256 × 256	2.0732E − 3	1.9839	1.4017	0.9822
Incomplete PIFE (ε = 0)	8 × 8	1.0751	-	23.1839	-
	16 × 16	3.9068E − 1	1.4603	16.5401	0.4872
	32 × 32	1.1601E − 1	1.7518	9.7531	0.7620
	64 × 64	3.1393E − 2	1.8857	5.2968	0.8808
	128 × 128	8.2042E − 3	1.9360	2.7677	0.9364
	256 × 256	2.0726E − 3	1.9849	1.4011	0.9821
Symmetric PIFE (ε = −1)	8 × 8	1.0738	-	23.1231	-
	16 × 16	3.9080E − 1	1.4582	16.5280	0.4844
	32 × 32	1.1596E − 1	1.7528	9.7468	0.7619
	64 × 64	3.1392E − 2	1.8852	5.2953	0.8802
	128 × 128	8.2055E − 3	1.9357	2.7674	0.9362
	256 × 256	2.0730E − 3	1.9848	1.4011	0.9820

Finally, we would like to point out that, even for a larger jump of the diffusion coefficient, our PIFE method is also valid and possesses optimal-order convergence properties. We present numerical results for , in Table 7 as an example for illustration.

Table 7

Numerical results for ,

Scheme	Grid	∥u−u_h∥_0,Ω	Order	|u−u_h|_1,h	Order
Nonsymmetric PIFE (ε = 1)	8 × 8	1.2074E − 1	-	2.3185	-
	16 × 16	3.1092E − 2	1.9573	1.1779	0.9769
	32 × 32	7.8627E − 3	1.9835	5.9427E − 1	0.9871
	64 × 64	1.9738E − 3	1.9940	2.9834E − 1	0.9942
	128 × 128	4.9306E − 4	2.0011	1.4939E − 1	0.9979
	256 × 256	1.2335E − 4	1.9991	7.4748E − 2	0.9990
Incomplete PIFE (ε = 0)	8 × 8	1.1978E − 1	-	2.3039	-
	16 × 16	3.1017E − 2	1.9492	1.1766	0.9695
	32 × 32	7.8636E − 3	1.9798	5.9412E − 1	0.9858
	64 × 64	1.9749E − 3	1.9934	2.9833E − 1	0.9939
	128 × 128	4.9317E − 4	2.0016	1.4936E − 1	0.9981
	256 × 256	1.2339E − 4	1.9989	7.4745E − 2	0.9987
Symmetric PIFE (ε = −1)	8 × 8	1.1967E − 1	-	2.3025	-
	16 × 16	3.1010E − 2	1.9482	1.1765	0.9687
	32 × 32	7.8665E − 3	1.9790	5.9412E − 1	0.9857
	64 × 64	1.9752E − 3	1.9937	2.9833E − 1	0.9938
	128 × 128	4.9317E − 4	2.0018	1.4936E − 1	0.9981
	256 × 256	1.2337E − 4	1.9991	7.4745E − 2	0.9987

Conclusions

In this paper, on non-body fitted triangular meshes, we have developed a partially penalty immersed Crouzeix-Raviart finite element method for both isotropic and anisotropic elliptic interface problems. For the theoretical analysis, we restrict ourselves to the isotropic interface problems for simplicity. The solvability and optimal error estimates are obtained. Then numerical experiments for both isotropic and anisotropic problems are conducted to illustrate the validity of our schemes.

In this method, edge averages on three edges are used as degrees of freedom, which makes it more advantageous for solving anisotropic interface problems than the use of values at three vertices. From numerical examples in the above section, we can observe that it works well for the anisotropic problems, even if the entries of the tensor-coefficients () are negative, while in [23], by three counter-examples, it is pointed out that, for some specially selected entries of the definite-positive diffusion matrix and the intersection points of the interface with the edges on the interface elements, the piecewise linear Lagrange-nodal-polynomial satisfying the jump conditions cannot be uniquely determined by its values at three vertices.