1 Introduction

Physical phenomena in a domain consisting of multiple materials and/or multiphase fluids, which are immiscible and are divided by distinct interfaces, are often modeled by either identical or different partial differential equations with discontinuous coefficients on both sides of interfaces. These problems are generally called interface problems, sometimes called interaction problems in some specific scenarios such as fluid-structure interaction (FSI) problems, e.g., see [5, 17] and others references therein. In the past several decades, two major numerical approaches – the body-fitted mesh method and the body-unfitted mesh method – have been developed for tackling interface problems, which are classified by how the computational mesh and then the interface conditions are handled along the interface. In contrast to the body-fitted mesh method such as the arbitrary Lagrangian-Eulerian (ALE) method [8, 9], which is obliged to adapt the mesh to accommodate the motion of the interface, the body-unfitted mesh method, due to its simplicity in the mesh generation, becomes more promising and more advantageous for interface problems whose interfaces may bear a large deformation/displacement, such as the immersed boundary method (IBM) [12, 14], the distribute Lagrange multiplier/fictitious domain (DLM/FD) method [1, 7], the immersed finite element method (IFEM) [10, 11], and etc.

Among the aforementioned body-unfitted mesh methods, taking all properties of reliability, accuracy, flexibility and theoretical guarantees into consideration, the DLM/FD method has shown a lot of strengths and potentials in theoretical analyses as well as practical applications for general interface problems, and has gained considerable popularity in simulating FSI problems as well. So in this paper, we focus on the DLM/FD method, where a fictitious equation that is defined in one subdomain is introduced to cover the other subdomain, and its mesh is fixed in the entire domain as a background mesh and needs not to be updated even if the interface moves or deforms. Benefited from this feature, the DLM/FD method has became more popular in the simulation of FSI problems, especially in the case of an immersed structure with large deformation/displacements. To enforce the interface conditions, the DLM/FD method introduces the Lagrange multiplier (a pseudo body force) to weakly enforce the fictitious variables equal everywhere to the primary variables of the equation defined in the immersed domain and on the interface too. A monolithic system bearing a saddle-point structure is thus formed in regard to the Lagrange multiplier and primary variables. Therefore, the classical Babus̆ka–Brezzi’s theory [3, 4] can be employed to prove the well-posedness, stability as well as convergence properties of the DLM/FD finite element method [2, 15].

To implement the DLM/FD method, an accurate and also efficient numerical integration scheme is needed to calculate the integration in which the distributed Lagrange multiplier is involved. For instance, (21) can be referred to in advance to preview the significance, where, in the finite element computation on each immersed element for the Lagrange multiplier terms denoted by the dual inner product \(\left\langle \lambda ,v_{h}|_{\varOmega }\right\rangle _{\varOmega _2}\), the integrand function is a product of two piecewise polynomials which are defined on two non-matched meshes, \(T_h(\varOmega )\) and \(T_H(\varOmega _2)\). Although the piecewise polynomial defined on \(T_h(\varOmega )\), \(v_h|_{\varOmega }\), can be transferred to \(T_H(\varOmega _2)\) through the interpolation approach, it is no longer sufficiently smooth in each immersed element of \(T_H(\varOmega _2)\) just because of the non-matching between \(T_h(\varOmega )\) and \(T_H(\varOmega _2)\). Thus we can not conclude that the commonly-used higher order numerical integration scheme will lead to a higher accuracy for those Lagrange multiplier terms.

The object of this paper is to study the effects of various numerical integration schemes on the performance of the DLM/FD methods for solving different interface problems with jump coefficients. Essentially, we do not want to let the accuracy of numerical integration influence the overall approximation accuracy, especially when the DLM/FD method is used for simulating complex problems, e.g., FSI problems. However, we find out different numerical integration methods indeed have significant effects on approximation errors of the DLM/FD method for different types of interface problems. Three types of numerical integration schemes are considered in our study: the commonly-used numerical integration schemes (see e.g. [6]), the compound type formulas, and the subgrid integration scheme proposed in [18]. Numerical results presented in Sect. 4 shows that the performance of DLM/FD method for solving elliptic interface problems is insensitive to the numerical integration but is sensitive to Stokes- and Stokes/elliptic- interface problems, and that the subgrid integration scheme always leads to the approximation solution of the highest accuracy.

The rest of this paper is organized as follows. The DLM/FD finite element method for solving different types of interface problems are recalled in Sect. 2. Several commonly-used numerical integration schemes, the compound type formulas and the subgrid integration technique are introduced in Sect. 3. Numerical performances are shown in Sect. 4.

2 DLM/FD Method for Three Types of Interface Problems

In this section, we briefly recall the DLM/FD finite element method for solving three different types of interface problems with jump coefficients.

Fig. 1.
figure 1

Graphical depiction of the domain with an immersed interface.

Full size image

2.1 Three Types of Interface Problems

The first type is the elliptic interface problem, defined as

$$\begin{aligned} -\nabla \cdot (\beta _1\nabla u_1)&= f_1,\qquad \quad \text {in}\ \varOmega _1,\end{aligned}$$
(1)
$$\begin{aligned} -\nabla \cdot (\beta _2\nabla u_2)&= f_2,\qquad \quad \text {in}\ \varOmega _2,\end{aligned}$$
(2)
$$\begin{aligned} u_1&=u_2,\qquad \quad \text {on}\ \varGamma ,\end{aligned}$$
(3)
$$\begin{aligned} \beta _1\nabla u_1\cdot \varvec{n}_1+\beta _2\nabla u_2\cdot \varvec{n}_2&= w,\qquad \quad \; \text {on}\ \varGamma ,\end{aligned}$$
(4)
$$\begin{aligned} u_1=0,\, u_2&=0,\qquad \quad \; \text {on}\ \partial \varOmega \backslash \varGamma , \end{aligned}$$
(5)

where, \(\varOmega = \varOmega _1\cup \varOmega _2 \subset \mathcal {R}^d\) (see Fig. 1), the interface \(\varGamma =\partial \varOmega _2\) is a closed curve that divides the domain \(\varOmega \) into an interior region \(\varOmega _2\) and an exterior region \(\varOmega _1\), \(\varvec{n}_1\) and \(\varvec{n}_2\) stand for the unit outward normal vectors on \(\partial \varOmega _1\) and \(\partial \varOmega _2\), respectively. u, that is defined in \(\varOmega \), satisfies \(u|_{\varOmega _1}=u_1,\ u|_{\varOmega _2}=u_2\) which are associated with \(f\in L^2(\varOmega )\) and \(f|_{\varOmega _1}=f_1\in L^2(\varOmega _1)\), \(f|_{\varOmega _2}=f_2\in L^2(\varOmega _2)\). \(\beta \in L^\infty (\varOmega )\) satisfies \(\beta |_{\varOmega _1}=\beta _1\in W^{1,\infty }(\varOmega _1),\ \beta |_{\varOmega _2}=\beta _2\in W^{1,\infty } (\varOmega _2)\) and \(\beta _1\ne \beta _2\).

The second type is the Stokes interface problem, defined as

$$\begin{aligned} -\nabla \cdot (\beta _1\nabla \varvec{u}_1)+\nabla p_1&= \varvec{f}_1,\qquad \quad \; \text {in}\ \varOmega _1,\end{aligned}$$
(6)
$$\begin{aligned} \nabla \cdot \varvec{u}_1&=0,\qquad \quad \;\;\; \text {in}\ \varOmega _1,\end{aligned}$$
(7)
$$\begin{aligned} -\nabla \cdot (\beta _2\nabla \varvec{u}_2)+\nabla p_2&= \varvec{f}_2,\qquad \quad \; \text {in}\ \varOmega _2,\end{aligned}$$
(8)
$$\begin{aligned} \nabla \cdot \varvec{u}_2&=0,\qquad \quad \;\;\; \text {in}\ \varOmega _2,\end{aligned}$$
(9)
$$\begin{aligned} \varvec{u}_1&=\varvec{u}_2,\qquad \quad \; \text {on}\ \varGamma ,\end{aligned}$$
(10)
$$\begin{aligned} (\beta _1\nabla \varvec{u}_1-p_1\varvec{I})\varvec{n}_1+(\beta _2\nabla \varvec{u}_2-p_2\varvec{I})\varvec{n}_2&= \varvec{w},\qquad \quad \;\; \text {on}\ \varGamma ,\end{aligned}$$
(11)
$$\begin{aligned} \varvec{u}_1=0,\, \varvec{u}_2&=0,\qquad \quad \;\;\; \text {on}\ \partial \varOmega \backslash \varGamma , \end{aligned}$$
(12)

which can be considered as the linearized model of an immiscible two-phase fluid flow problem, where the primary variable pair (\(\varvec{u},p\)) satisfies \(\varvec{u}|_{\varOmega _1}=\varvec{u}_1,\ \varvec{u}|_{\varOmega _2}=\varvec{u}_2,\ p|_{\varOmega _1}=p_1,\ p|_{\varOmega _2}=p_2\), which are associated with \(\varvec{f}\in (L^2(\varOmega ))^d\) satisfying \(\varvec{f}|_{\varOmega _1}=\varvec{f}_1\in (L^2(\varOmega _1))^d\), \(\varvec{f}|_{\varOmega _2}=\varvec{f}_2\in (L^2(\varOmega _2))^d\). The jump coefficient \(\beta \in L^\infty (\varOmega )\) satisfies \(\beta |_{\varOmega _1}=\beta _1\in W^{1,\infty } (\varOmega _1),\ \beta |_{\varOmega _2}=\beta _2\in W^{1,\infty } (\varOmega _2),\ \beta _1\ne \beta _2\).

And, the third type is the Stokes/elliptic interface problem, defined as

$$\begin{aligned} -\nabla \cdot (\beta _1\nabla \varvec{u}_1)+\nabla p_1&= \varvec{f}_1,\qquad \quad \; \text {in}\ \varOmega _1,\end{aligned}$$
(13)
$$\begin{aligned} \nabla \cdot \varvec{u}_1&=0,\qquad \quad \;\;\; \text {in}\ \varOmega _1,\end{aligned}$$
(14)
$$\begin{aligned} -\nabla \cdot (\beta _2\nabla \varvec{u}_2)&= \varvec{f}_2,\qquad \quad \; \text {in}\ \varOmega _2,\end{aligned}$$
(15)
$$\begin{aligned} \varvec{u}_1&=\varvec{u}_2,\qquad \quad \; \text {on}\ \varGamma ,\end{aligned}$$
(16)
$$\begin{aligned} (\beta _1\nabla \varvec{u}_1-p_1\varvec{I})\varvec{n}_1+\beta _2\nabla \varvec{u}_2\varvec{n}_2&= \varvec{w},\qquad \quad \;\; \text {on}\ \varGamma ,\end{aligned}$$
(17)
$$\begin{aligned} \varvec{u}_1=0,\, \varvec{u}_2&=0,\qquad \quad \;\;\; \text {on}\ \partial \varOmega \backslash \varGamma , \end{aligned}$$
(18)

which can be considered as a steady state linearized model of FSI problems, where, the primary variable pair (\(\varvec{u},p_1\)) satisfies \(\varvec{u}|_{\varOmega _1}=\varvec{u}_1,\ \varvec{u}|_{\varOmega _2}=\varvec{u}_2,\ p_1\in \varOmega _1\), which are associated with \(\varvec{f}\in (L^2(\varOmega ))^d\) satisfying \(\varvec{f}|_{\varOmega _1}=\varvec{f}_1\in (L^2(\varOmega _1))^d\), \(\varvec{f}|_{\varOmega _2}=\varvec{f}_2\in (L^2(\varOmega _2))^d\). The jump coefficient \(\beta \in L^\infty (\varOmega )\) satisfies \(\beta |_{\varOmega _1}=\beta _1\in W^{1,\infty } (\varOmega _1),\ \beta |_{\varOmega _2}=\beta _2\in W^{1,\infty } (\varOmega _2),\ \beta _1\ne \beta _2\).

2.2 DLM/FD Method for the Elliptic Interface Problem

The main idea of the DLM/FD method is to smoothly extend one material, such as the fluid in FSI, into another material’s subdomain, such as the structure in FSI, as a fictitious equation whose primary variable is constrained to equal to that of the occupied material’s equation in its subdomain and on the interface too. And, such constraint is weakly imposed through the Lagrange multiplier (a pseudo body force) in the DLM/FD method. Thus, a monolithic saddle-point system is formed in regard to the Lagrange multiplier and primary variables, for which the classical Babus̆ka–Brezzi’s theory [3, 4] can be employed to prove the well-posedness, stability as well as convergence properties [2, 16]. In the following, we introduce the DLM/FD finite element method for each type of interface problem that is defined in Sect. 2.1.

Define \(\varvec{V}^E =H^1_0(\varOmega ),\quad \varvec{V}^E_2 = H^1(\varOmega _2),\quad \varvec{\varLambda }^E = \left( \varvec{V}^E_2\right) ^*,\) where \(\left( \varvec{V}^E_2\right) ^*\) denotes the dual space of \(\varvec{V}^E_2\). Let \(T_h(\varOmega )\) and \(T_H(\varOmega _2)\) be the meshes of \(\varOmega \) and \(\varOmega _2\), respectively. And denote by \(\varvec{V}^E_h\), \(\varvec{V}^E_{2,H}\) and \(\varvec{\varLambda }^E_H\) the conforming finite element spaces of \(\varvec{V}^E\), \(\varvec{V}^E_2\) and \( {\varvec{\varLambda }^E}\), respectively. Here and hereafter, let \((\cdot ,\cdot )_{\omega }\) be the \(L^2\) inner product over \(\omega \), and \(\left\langle \cdot ,\cdot \right\rangle _{\omega }\) be the dual product over \(\omega \) that is actually associated with the \(H^1\) inner product [1].

Then, the DLM/FD finite element method for solving elliptic interface problem (1)–(5) can be defined as follows [1, 2]: Find \((\tilde{u}_h,\ u_{2,H}, \ \lambda _H)\in \varvec{V}^E_h\times \varvec{V}^E_{2,H} \times {\varvec{\varLambda }^E_H}\) such that

$$\begin{aligned} (\tilde{\beta }\nabla \tilde{u}_{h},\nabla v_{h})_\varOmega +\left\langle \lambda _{H},v_{h}|_{\varOmega _2}\right\rangle _{\varOmega _2}=(\tilde{f}, v_{h})_\varOmega , \end{aligned}$$
(19)
$$\begin{aligned} \left( (\beta _2-\tilde{\beta })\nabla u_{2,H},\nabla v_{2,H}\right) _{\varOmega _2}-\left\langle \lambda _{H}, v_{2,H}\right\rangle _{\varOmega _2} =\left( f_2-\tilde{f}_{2}, v_{2,H}\right) _{\varOmega _{2}}+(w,v_{2,H})_{\varGamma }, \end{aligned}$$
(20)
$$\begin{aligned} \left\langle \xi _{H},\tilde{u}_{h}|_{\varOmega _2}- u_{2,H}\right\rangle _{\varOmega _2}=0, \forall ~ ( v_{h}, \ v_{2,H},\ \xi ) \in \varvec{V}^{E}_{h}\times \varvec{V}^{E}_{2,H}\times {\varvec{\varLambda }^E_H}. \end{aligned}$$
(21)

2.3 DLM/FD Method for the Stokes Interface Problem

Define \(\varvec{V}^S =\left( H^1_0(\varOmega )\right) ^d,\quad \varvec{Q}^S = L^2(\varOmega ),\quad \varvec{V}^S_2 = \left( H^1(\varOmega _2)\right) ^d,\quad \varvec{\varLambda }^S = \left( \varvec{V}^S_2\right) ^*\), where \(\left( \varvec{V}^S_2\right) ^*\) denotes the dual space of \(\varvec{V}^S_2\). Let \(\varvec{V}^S_h\), \(\varvec{V}^S_{2,H}\), \(\varvec{Q}^S_h\) and \(\varvec{\varLambda }^S_H\) be the conforming finite element spaces of \(\varvec{V}^S\), \(\varvec{V}^S_2\), \(\varvec{Q}^S\) and \( {\varvec{\varLambda }^S}\), respectively. Then, the DLM/FD finite element method for solving the Stokes interface problem (6)–(12) can be defined as follows [13]: Find \((\tilde{\varvec{u}}_h,\ \varvec{u}_{2,H}, \ \tilde{p}_h ,\ \varvec{\lambda }_H)\in \varvec{V}^S_h\times \varvec{V}^S_{2,H}\times Q^S_h\times {\varvec{\varLambda }^S_H}\) such that

$$\begin{aligned} (\tilde{\beta }\nabla \tilde{\varvec{u}}_h,\nabla \varvec{v}_h)_\varOmega -(\tilde{p}_h,\nabla \cdot \varvec{v}_h)_\varOmega +\left\langle \varvec{\lambda }_H,\varvec{v}_h|_{\varOmega _2}\right\rangle _{\varOmega _2}=(\tilde{\varvec{f}},\varvec{v}_h)_\varOmega , \end{aligned}$$
(22)
$$\begin{aligned} (\nabla \cdot \tilde{\varvec{u}}_h,q_h)_\varOmega =0,\end{aligned}$$
(23)
$$\begin{aligned} \left( (\beta _2-\tilde{\beta })\nabla \varvec{u}_{2,H},\nabla \varvec{v}_{2,H}\right) _{\varOmega _2}-\left\langle \varvec{\lambda }_H,\varvec{v}_{2,H}\right\rangle _{\varOmega _2} =\left( \varvec{f}_2-\tilde{\varvec{f}}_2,\varvec{v}_{2,H}\right) _{\varOmega _2}+(\varvec{w},\varvec{v}_{2,H})_{\varGamma }, \end{aligned}$$
(24)
$$\begin{aligned} \left\langle \varvec{\xi }_H,\tilde{\varvec{u}}_h|_{\varOmega _2}-\varvec{u}_{2,H}\right\rangle _{\varOmega _2}=0, \forall \ (\varvec{v}_h, \ \varvec{v}_{2,H},\ q_h,\ \varvec{\xi }_H) \in \varvec{V}^S_h\times \varvec{V}_{2,H}^S\times Q^S_h\times {\varvec{\varLambda }^S_H}. \end{aligned}$$
(25)

2.4 DLM/FD Method for the Stokes/Elliptic Interface Problem

The spaces to be used for defining the weak formulation and the DLM/FD formulation of Stokes/elliptic interface problem are the same as those of Stokes interface problem. The DLM/FD finite element method is proposed and analyzed in [15] for the Stokes/elliptic interface problem (13)–(18), as described below.

Find \((\tilde{\varvec{u}}_h,\ \varvec{u}_{2,H}, \ \tilde{p}_h,\ \varvec{\lambda }_H)\in \varvec{V}^S_h\times \varvec{V}^S_{2,H}\times Q^S_h\times {\varvec{\varLambda }^S_H}\) such that

$$\begin{aligned} (\tilde{\beta }\nabla \tilde{\varvec{u}}_h,\nabla \varvec{v}_h)_\varOmega -(\tilde{p}_h,\nabla \cdot \varvec{v}_h)_\varOmega +\left\langle \varvec{\lambda }_H,\varvec{v}_h|_{\varOmega _2}\right\rangle _{\varOmega _2}=(\tilde{\varvec{f}},\varvec{v}_h)_\varOmega ,\end{aligned}$$
(26)
$$\begin{aligned} (\nabla \cdot \tilde{\varvec{u}}_h,q_h)_\varOmega -(\nabla \cdot \varvec{u}_{2,H},q_h)_{\varOmega _2}=0,\end{aligned}$$
(27)
$$\begin{aligned} \left( (\beta _2-\tilde{\beta })\nabla \varvec{u}_{2,H},\nabla \varvec{v}_{2,H}\right) _{\varOmega _2}+\left( p_h|_{\varOmega _2},\nabla \cdot \varvec{v}_{2,H}\right) _{\varOmega _2}\nonumber \\ -\left\langle \varvec{\lambda }_H,\varvec{v}_{2,H}\right\rangle _{\varOmega _2} =\left( \varvec{f}_2-\tilde{\varvec{f}}|_{\varOmega _2},\varvec{v}_{2,H}\right) _{\varOmega _2}+(\varvec{w},\varvec{v}_{2,H})_{\varGamma },\end{aligned}$$
(28)
$$\begin{aligned} \left\langle \varvec{\xi }_H,\tilde{\varvec{u}}_h|_{\varOmega _2}-\varvec{u}_{2,H}\right\rangle _{\varOmega _2}=0, \forall \ (\varvec{v}_h,\ \varvec{v}_{2,H}, \ q_h,\ \varvec{\xi }_H)\in \varvec{V}^S_h\times \varvec{V}^S_{2,H}\times Q^S_h\times {\varvec{\varLambda }^S_H}. \end{aligned}$$
(29)

3 Numerical Integration Schemes

Note that in the DLM/FD finite element methods described in Sects. 2.2, 2.3 and 2.4, some inner product terms involve an integration over the immersed domain \(\varOmega _2\) with a product of two piecewise functions defined on \(T_h(\varOmega )\) and \(T_H(\varOmega _2)\), respectively. For instance,

$$\begin{aligned}&\quad \left\langle \lambda _H,v_h|_{\varOmega _2}\right\rangle _{\varOmega _2} \text { in } (19) \text { and } \left\langle \xi _H,\tilde{u}_h|_{\varOmega _2}\right\rangle _{\varOmega _2} \text { in } (21), \end{aligned}$$
(30)
$$\begin{aligned}&\quad \left\langle \varvec{\lambda }_H,\varvec{v}_h|_{\varOmega _2}\right\rangle _{\varOmega _2} \text { in } (22) \text { and } \left\langle \varvec{\xi }_H,\tilde{\varvec{u}}_h|_{\varOmega _2}\right\rangle _{\varOmega _2}\text { in } (25),\end{aligned}$$
(31)
$$\begin{aligned}&\quad \left( \nabla \cdot \tilde{\varvec{u}}_{2,H},q_h\right) _{\varOmega _2}\text { in } (26)\text { and } \left( p_h|_{\varOmega _2},\nabla \cdot \varvec{v}_{2,H}\right) _{\varOmega _2}\text { in } (29). \end{aligned}$$
(32)

Since two grids \(T_h(\varOmega )\) and \(T_H(\varOmega _2)\) are constructed independently, they are generally not matched with each other, see Fig. 2. Thus, we need to employ some numerical techniques to implement the numerical integration for terms shown in (30)–(32) in order to make the DLM/FD finite element method perform well. In this paper, we restrict ourselves to the cases of triangular finite element.

It shall be pointed out that we can not simply make such a conclusion that the higher order the integration scheme, the higher accuracy the numerical integration. That is because the integrand functions shown in (30)–(32) are essentially the product of piecewise polynomials defined on two non-matched meshes, \(T_h(\varOmega )\) or \(T_H(\varOmega _2)\), inducing an insufficiently smooth integrand function in each element of \(T_H(\varOmega _2)\), and thus it is not able to deduce a higher order derivative in the remainder of the numerical integration scheme.

3.1 Commonly-Used Numerical Integration Schemes

The numerical integration schemes in the field of finite element method are of the form: \(\int _{\hat{K}} f(\hat{x}) d\hat{x} = \sum _{i} \omega _i f(\lambda _i)\), where \(\hat{K}\) stands for the reference element that is an isosceles right triangle with two sides equal 1, that is, three vertices of \(\hat{K}\) are (0, 0), (1, 0) and (0, 1). And, the quadrature points are denoted by the barycentric coordinates. Nine different commonly-used numerical integration schemes, from Scheme 1 to Scheme 9 as listed in the Appendix, are adopted to carry out our numerical studies in Sect. 4.

Moreover, to improve the numerical integration accuracy, we consider the compound type formulas of numerical integration, for which each triangular reference element is divided into four equilateral sub-triangles. For example, Scheme 4c1 can be obtained by applying Scheme 4 of the numerical integration to each one of four sub-triangles in the first-time refinement of the reference element, Scheme 4c2 is constructed by applying Scheme 4 to each one of 16 sub-triangles after refining the reference element twice, and Scheme 4c3 based on Scheme 4 and 64 sub-triangles after refining the reference element for three times, and so forth, which are not necessarily shown in the Appendix for the simplicity.

3.2 The Subgrid Integration Scheme

The subgrid integration scheme has been proposed and used in [18] for solving two-dimensional parabolic interface problems with the DLM/FD method. The main ingredient of this scheme is to generate a subgrid, \(T^r_H(\varOmega _2)\), by finding all intersections of \(T_h(\varOmega )\) and \(T_H(\varOmega _2)\) then forming a new subgrid structure in the immersed domain \(\varOmega _2\). Clearly, \(T^r_H(\varOmega _2)\) is a matching subset of both \(T_h(\varOmega )\) and \(T_H(\varOmega _2)\), and \(T^r_H(\varOmega _2)\) is a finer mesh on \(\varOmega _2\). Thus, all integral terms as shown in (30)–(32) can be implemented on \(T^r_H(\varOmega _2)\) now, and the piecewise polynomials defined on \(T_h(\varOmega )\) and \(T_H(\varOmega _2)\) must be smooth in each element of \(T^r_H(\varOmega _2)\). Then, the interpolations between two different meshes are avoided. An example of such a subgrid is shown in Fig. 2.

Since the integrand function on each element of the subgrid \(T^r_H(\varOmega _2)\) is smooth, the commonly-used integration schemes with the help of such a subgrid can yield sufficiently high accuracy for numerical integrations. In this sense, the subgrid integration scheme is of the highest accuracy among the schemes presented in this paper, as long as a subtle implementation programming can be done to eliminate every possible geometrical error when finding the intersection points of two non-matching meshes. Interested readers can refer to [18] for more details about the subgrid integration scheme and its implementation.

4 Numerical Experiments

In what follows, we will carry out two scenarios in our numerical experiments to investigate effects of different numerical integration schemes on the approximation solutions of the DLM/FD finite element method for solving three types of steady interface problems: (1) the real solution of the original interface problem is unknown; (2) the real solution of the original interface problem is known.

To that end, a fixed mesh size is chosen as \(h=1/64\) and \(H=1/64\) to attain two meshes \(T_h(\varOmega )\) and \(T_H(\varOmega _2)\), as depicted in Fig. 2. On these two fixed meshes, we will implement the DLM/FD finite element method with different numerical integration schemes for the elliptic-, Stokes-, and Stokes/elliptic-interface problems, respectively.

Fig. 2.
figure 2

The meshes (from left to right): the background mesh \(T_h(\varOmega )\) partitioned in a square \(\varOmega \), the foreground mesh \(T_H(\varOmega _2)\) partitioned in a disk \(\varOmega _2\), and the subgrid \(T^r_H(\varOmega _2)\) in \(\varOmega _2\).

Full size image

4.1 The Scenario of Unknown Real Solutions

In this scenario, since the real solution is unknown, we will report the energy norm of the numerical solution instead of its approximation error to the real solution, and compare all numerical solutions in the energy norm obtained from different numerical integration schemes for each type of interface problem.

The Elliptic Interface Problem with a Unknown Real Solution. Let \(\varOmega \) be a square with the size \((0,6)\times (0,6)\) and \(\varOmega _2\) be a unit circle with the center at (3, 3). Consider the elliptic interface problem (1)–(5) with the following jump coefficients and right hand side functions:

  • Case E1: \(\beta _1 = 1\), \(\beta _2=100\), \(f_1=f_2=1\),

  • Case E2: \(\beta _1 = 1\), \(\beta _2=10000\), \(f_1=f_2=1\).

The linear (\(P_1\)) finite element space is adopted to define \(\varvec{V}^E_h\), \(\varvec{V}^E_{2,H}\) and \(\varvec{\varLambda }^E_H\), which are used in the implementation of the DLM/FD finite element method (19)–(21). Numerical solutions in \(H^1\)-norm in regard to different numerical integration schemes are reported in Tables 1 and 2 for two cases E1–E2, where, we use the numbers from 1 to 9 to successively represent nine commonly-used non-compound type numerical integration schemes, use the notation “4c1”, “4c2” and “4c3” to represent three compound-type formula of numerical integration, respectively, and use “s” to stand for the subgrid integration scheme.

Table 1. Numerical results of Case E1
Full size table
Table 2. Numerical results of Case E2
Full size table

From Tables 1 and 2, we find out all numerical integration schemes yield almost the same numerical solution in terms of the energy (\(H^1\)-) norm, that is even true for the case of large jump coefficients. Thus, the DLM/FD finite element method for the elliptic interface problem is insensitive with different numerical integration schemes, i.e., all commonly-used numerical integration schemes and their compound formulas can achieve the same approximation accuracy with the subgrid integration scheme which supposes to produce the highest accuracy of the numerical integration for the DLM/FD method.

The Stokes Interface Problem with a Unknown Real Solution. Let \(\varOmega \) be a square with the size \((0,1)\times (0,1)\) and \(\varOmega _2\) be a circle with the center at (0.3, 0.3) and the radius 0.1. Consider the Stokes interface problem (6)–(12) with the following jump coefficients and right hand side functions:

  • Case S1: \(\beta _1 = 1\), \(\beta _2=100\), \(\varvec{f}_1=(e^x,e^y)^T\), \(\varvec{f}_2=(0,0)^T\).

  • Case S2: \(\beta _1 = 1\), \(\beta _2=10000\), \(\varvec{f}_1=(e^x,e^y)^T\), \(\varvec{f}_2=(0,0)^T\).

Table 3. Numerical results of Case S1
Full size table
Table 4. Numerical results of Case S2
Full size table

The Taylor-Hood \(P_2P_1\)-mixed finite element is adopted to define \(\varvec{V}^S_h\), \(\varvec{V}^S_{2,H}\), \(\varvec{Q}^S_h\) and \(\varvec{\varLambda }^S_H\), which are used in the implementation of the DLM/FD finite element method (22)–(25). Numerical solutions in the energy norm (the velocity in \(H^1\)-norm and the pressure in \(L^2\)-norm) in regard to different numerical integration schemes are reported in Tables 3 and 4 for two cases S1–S2. From these tables, we observe some obvious differences on the numerical solutions of velocity between different numerical integration schemes, while the numerical pressure is insensitive with different integration schemes.

The Stokes/Elliptic Interface Problem with a Unknown Real Solution. Let \(\varOmega \) be a square with the size \((0,1)\times (0,1)\) and \(\varOmega _2\) be a circle with the center at (0.3, 0.3) and the radius 0.1. Consider the Stokes/elliptic interface problem (13)–(18) with the following jump coefficients and right hand side functions:

  • Case SE1: \(\beta _1 = 1\), \(\beta _2=100\), \(\varvec{f}_1=\varvec{f}_2=(1,1)^T\),

  • Case SE2: \(\beta _1 = 1\), \(\beta _2=10000\), \(\varvec{f}_1=\varvec{f}_2=(1,1)^T\).

The same \(P_2P_1\)-mixed finite element space is adopted to define \(\varvec{V}^S_h\), \(\varvec{V}^S_{2,H}\), \(\varvec{Q}^S_h\) and \(\varvec{\varLambda }^S_H\), which are used in the implementation of the DLM/FD finite element method (26)–(29). Numerical solutions in the energy norm in regard to different numerical integration schemes are reported in Tables 5 and 6 for two cases SE1–SE2, from which the similar conclusions can be drawn as those for the numerical results of Stokes interface problem in Cases S1–S2.

Table 5. Numerical results of Case SE1
Full size table
Table 6. Numerical results of Case SE2
Full size table

4.2 The Scenario of Known and Smooth Real Solutions

The Elliptic Interface Problem with a Known Real Solution. Let \(\varOmega \) be a square with the size \((0,6)\times (0,6)\) and \(\varOmega _2\) be a unit circle with the center at (3, 3). Consider the elliptic interface problem (1)–(5) with the following jump coefficients and right hand side functions:

  • Case E3: \(\beta _1 = 1\), \(\beta _2=100\), \(f_1\) and \(f_2\) are defined below;

  • Case E4: \(\beta _1 = 1\), \(\beta _2=10000\), \(f_1\) and \(f_2\) are defined below,

where, \(f_1\) and \(f_2\) are chosen such that the real solution of (1)–(5) is taken as \(u = (x-6)x(y-6)y((x-3)^2+(y-3)^2-1)^2.\)

Table 7. Numerical results of Case E3
Full size table
Table 8. Numerical results of Case E4
Full size table

With the \(P_1\) finite element space, the numerical results of the DLM/FD finite element method (19)–(21) are reported in Tables 7 and 8 for two cases E3-E4, from which we observe that all numerical errors in \(H^1\)-norm in regard to different numerical integration schemes are almost the same, confirming again that the DLM/FD finite element method for the elliptic interface problem is insensitive with different numerical integration schemes.

The Stokes Interface Problem with a Known Real Solution. Let \(\varOmega \) be a square with the size \((0,1)\times (0,1)\) and \(\varOmega _2\) be a circle with the center at (0.3, 0.3) and the radius 0.1. Consider the Stokes interface problem (6)–(12) with the following jump coefficients and right hand side functions:

  • Case S3: \(\beta _1 = 1\), \(\beta _2=100\), \(\varvec{f}_1\) and \(\varvec{f}_2\) are defined below;

  • Case S4: \(\beta _1 = 1\), \(\beta _2=10000\), \(\varvec{f}_1\) and \(\varvec{f}_2\) are defined below,

where, \(\varvec{f}_1\) and \(\varvec{f}_2\) are chosen such that the real solution of (6)–(12) can be taken as

$$\begin{aligned} \varvec{u}&=\left( \begin{array}{l} \frac{(y-0.3)((x-0.3)^2+(y-0.3)^2-0.01)}{\beta _i}\\ \frac{-(x-0.3)((x-0.3)^2+(y-0.3)^2-0.01)}{\beta _i} \end{array} \right) \quad \text { if } (x,y)^T\in \varOmega _i,\\ p&= 0.01(x^3-y^3). \end{aligned}$$

Note that an inhomogeneous boundary condition is employed instead for the Stokes interface problem in the above two cases.

With the \(P_2P_1\)-mixed finite element space, the numerical results of the DLM/FD finite element method (22)–(25) are reported in Tables 9 and 10 for two cases S3-S4, from which we observe the obvious differences on numerical solutions of both the velocity and the pressure between different numerical integration schemes, and, the subgrid integration scheme produces the lowest approximation error, i.e., the best accuracy, for both the velocity and the pressure in comparison with other numerical integration schemes in all two cases.

To look into effects of other numerical integration schemes than the subgrid scheme, in terms of the following formula

$$\begin{aligned} \frac{\left| \Vert \varvec{u}-\tilde{\varvec{u}}_h^{i}\Vert _1-\Vert \varvec{u}-\tilde{\varvec{u}}_h^{s}\Vert _1\right| }{\Vert \varvec{u}-\tilde{\varvec{u}}_h^{s}\Vert _1} \text { or } \frac{\left| \Vert \varvec{u}-\tilde{\varvec{u}}_h^{5ci}\Vert _1-\Vert \varvec{u}-\tilde{\varvec{u}}_h^{s}\Vert _1\right| }{\Vert \varvec{u}-\tilde{\varvec{u}}_h^{s}\Vert _1},\\ \frac{\left| \Vert p-p_h^{i}\Vert _0-\Vert p-p_h^{s}\Vert _0\right| }{\Vert p-p_h^{s}\Vert _0} \text { or } \frac{\left| \Vert p-p_h^{5ci}\Vert _0-\Vert p-p_h^{s}\Vert _0\right| }{\Vert p-p_h^{s}\Vert _0}, \end{aligned}$$

we calculate and plot the relative differences between numerical errors obtained from other integration schemes and those from the subgrid scheme, as shown in Fig. 3, where, two cases illustrate a slow convergence tendency on numerical errors of the velocity towards those from the subgrid scheme along with the increase of numerical integration scheme number. However, such tendency does not even apply to numerical errors of the pressure, instead, Fig. 3 shows that a higher order integration scheme may lead to a worse approximation to the pressure. Furthermore, in Fig. 4, a special case is observed that the numerical solution obtained by Scheme 2 is significantly polluted inside \(\varOmega _2\) while the numerical solution obtained by the subgrid integration scheme is of better accuracy in \(\varOmega _2\), which implies that the Lagrange multiplier does not work well to enforce (25) hold true, reflecting from the low accuracy of Scheme 2.

Table 9. Numerical results of Case S3
Full size table
Table 10. Numerical results of Case S4
Full size table
Fig. 3.
figure 3

Relative differences of the velocity error in \(H^1\)-norm and of the pressure error in \(L^2\)-norm for Cases S3 and S4 (from left to right).

Full size image
Fig. 4.
figure 4

The first components of \(\varvec{u}-\tilde{\varvec{u}}^2_h\) and \(\varvec{u}_2-\varvec{u}^2_{2,H}\) vs the first components of \(\varvec{u}-\tilde{\varvec{u}}^s_h\) and \(\varvec{u}_2-\varvec{u}^2_{2,H}\) for Cases S3 (from left to right).

Full size image

The Stokes/Elliptic Interface Problem with a Known Real Solution. Let \(\varOmega \) be a square with the size \((0,1)\times (0,1)\) and \(\varOmega _2\) be a circle with the center at (0.3, 0.3) and the radius 0.1. Consider the Stokes/elliptic interface problem (13)–(18) with the following jump coefficients and right hand side functions:

  • Case SE3: \(\beta _1 = 1\), \(\beta _2=100\), \(\varvec{f}_1\) and \(\varvec{f}_2\) are defined below;

  • Case SE4: \(\beta _1 = 1\), \(\beta _2=10000\), \(\varvec{f}_1\) and \(\varvec{f}_2\) are defined below,

where, \(\varvec{f}_1\) and \(\varvec{f}_2\) are chosen such that the real solution of (13)–(18) is taken as

$$\begin{aligned} \varvec{u}&=\left( \begin{array}{l} \frac{(y-0.3)((x-0.3)^2+(y-0.3)^2-0.01)}{\beta _i}\\ \frac{-(x-0.3)((x-0.3)^2+(y-0.3)^2-0.01)}{\beta _i} \end{array} \right) \quad \text { if } (x,y)^T\in \varOmega _i,\\ p&= 0.01(x^3-y^3)((x-0.3)^2+(y-0.3)^2-0.01). \end{aligned}$$

Note that an inhomogeneous boundary condition is employed instead for the Stokes/elliptic interface problem in the above two cases.

With the \(P_2P_1\)-mixed finite element space, the numerical results of the DLM/FD finite element method (26)–(29) are reported in Tables 11 and 12 and Fig. 5 for two cases SE3–SE4, from which the same conclusions can be drawn as those for the numerical velocity errors of Stokes interface problem in Cases S3-S4, and, the slowly convergent relative differences of \(\Vert \tilde{\varvec{u}}-\tilde{\varvec{u}}_{h}\Vert _1\), that is generally between \(20\%\) and \(60\%\), can be observed in Fig. 5.

Table 11. Numerical results of Case SE3
Full size table
Table 12. Numerical results of Case SE4
Full size table

In summary, we notice that for the elliptic-, Stokes- and Stokes/elliptic- interface problem, their coefficient matrices of the linear algebra systems obtained by the DLM/FD FEM are in the following forms, respectively

$$\begin{aligned} S_E=\left( \begin{array}{ll|l} A&{}O&{}B^T\\ O&{}A_2&{}C^T\\ \hline B&{}C&{}O \end{array} \right) ,\ S_S=\left( \begin{array}{ll|ll} A&{}O&{}B^T&{}C^T\\ O&{}A_2&{}O&{}D^T\\ \hline B&{}O&{}O&{}O\\ C&{}D&{}O&{}O \end{array} \right) ,\ S_{SE}=\left( \begin{array}{ll|ll} A&{}O&{}B^T&{}C^T\\ O&{}A_2&{}E&{}D^T\\ \hline B&{}-E&{}O&{}O\\ C&{}D&{}O&{}O \end{array} \right) , \end{aligned}$$
(33)
Fig. 5.
figure 5

Relative differences of the velocity error in \(H^1\)-norm for Cases SE3 and SE4 (from left to right).

Full size image

where “O” denotes the zero matrix block. One way to view these linear algebra systems as saddle point problems is to split the matrices as showed in (33), where the Lagrange multiplier and the pressure are bundled in \(S_S\) and \(S_{SE}\). Slight changes of B in \(S_E\), C in \(S_S\), C and E in \(S_{SE}\) may result in quite different effects on the corresponding numerical solutions, which is a possible reason for explaining the observations illustrated in the previous sections: the DLM/FD FEM is insensitive with various numerical integration schemes for the elliptic interface problem, however, it is sensitive with different numerical integration schemes for both the Stokes- and the Stokes/elliptic interface problems. Although only two-dimensional interface problems are considered in the numerical experiments, the DLM/FD FEM and subgrid integration technique can be used for solving the interface problems in any dimension.

5 Conclusions

In this paper, we study the effects of different numerical integration schemes, including the commonly-used integration schemes and their compound-type schemes, and the subgrid integration scheme, on the performance of the DLM/FD finite element method for solving the elliptic-, Stokes-, and Stokes/elliptic-interface problem. Numerical experiments illustrate that: (1) DLM/FD FEM is insensitive with various numerical integration schemes for the elliptic interface problem; (2) DLM/FD FEM is sensitive with different numerical integration schemes for both Stokes- and Stokes/elliptic interface problem, and sometimes, numerical solutions obtained from other integration schemes than the subgrid integration scheme may not be reliable; (3) the subgrid integration scheme always results in numerical solutions with the best accuracy.