A NONCONFORMING IMMERSED VIRTUAL ELEMENT METHOD FOR ELLIPTIC INTERFACE PROBLEMS

. This paper presents the lowest-order nonconforming immersed virtual element method for solving elliptic interface problems on unfitted polygonal meshes. The local discrete space on each interface mesh element consists of the solutions of local interface problems with Neumann boundary conditions, and the elliptic projection is modified so that its range is the space of broken linear polynomials satisfying the interface conditions. We derive optimal error estimates in the broken 𝐻 1 -norm and 𝐿 2 -norm, under the piecewise 𝐻 2 -regulartiy assumption. In our scheme, the mesh assumptions for error analysis allow small cut elements. Several numerical experiments are provided to confirm the theoretical results


Introduction
In recent years, there has been a lot of interest in developing numerical methods for solving partial differential equations (PDEs) on general polygonal/polyhedral meshes, including mimetic finite difference (MFD) methods [10,19,20], hybrid high-order (HHO) methods [31,36], hybridizable discontinuous Galerkin (HDG) methods [32,33], weak Galerkin (WG) methods [60,65], and so on.Among them, the virtual element method (VEM), as an evolution of the MFD method into the framework of the finite element method (FEM), was introduced in [7].The main feature of the VEM is that the local shape functions, called the virtual elements, are defined implicitly as the solutions of certain local PDEs, and they are characterized by the degrees of freedom (DOFs).Although it is impossible to construct these functions explicitly in general, the VEM can be implemented using the DOFs only.The VEM also has been successfully developed for a wide range of problems: Stokes problem [11,24], elasticity problem [8,51,66], Maxwell problem [12,13], etc.We also refer to [2,5,9,21,49] and the references therein for more thorough survey.
On the other hand, there are numerous engineering and physical problems where the underlying PDEs have an interface, such as multiphase flows, solid mechanics with multiple materials, and Hele-Shaw flows, etc (see, e.g., [6,40,42,47]).The PDEs governing such problems involve with discontinuous coefficients across the interface, which usually leads to low global regularity of the solution, even when the interface is smooth.The low global regularity causes a deterioration in performance of the traditional FEMs, unless the mesh is aligned with the interface.However, it takes a lot of time to generate interface-fitted meshes when the interface is geometrically complicated or is moving as time evolves.For such cases, it may be more efficient to use unfitted or structured meshes than fitted meshes.Moreover, one can exploit geometric multigrid methods for the structured meshes.Researchers developed several numerical schemes using unfitted triangular or rectangular meshes: cut FEMs [39,40], extended FEMs [14,15], and immersed FEMs [52,55,57], to name just a few.In particular, the immersed FEM modifies the traditional finite elements so that they satisfy the interface conditions, while keeping the optimal approximation capabilities.Lagrange-type elements were studied in [41,55,57], while nonconformingtype elements were studied in [44,52,58].Other related works can be found in [27,38,43,45,48,50,53] and the references therein.
Due to the great flexibility of polygonal meshes in the mesh generation process, several researchers focused on developing interface-fitted polygonal mesh generators and analyzing schemes for interface problems on such meshes (see, e.g., [28,61,64]).Nevertheless, it would be still attractive to use unfitted polygonal meshes in some situations, such as problems involving moving interfaces as time evolves or during the computation of the free-boundary problems.Recently, several numerical schemes using unfitted polygonal meshes were developed.The authors in [22,23] proposed the unfitted HHO method for elliptic interface problems, where a Nitsche-type formulation is used.They proved that the method exhibits an optimal error estimate in the  1 -norm.However, to ensure the optimal convergence, it requires some additional mesh procedures, which prevent the appearance of small cut elements.The lowest-order Lagrange-type immersed VEM for triangular meshes was developed in [26].Unlike the Lagrange-type immersed FEM, the local shape functions are conforming, and DG-type consistency terms are not required to guarantee the optimal convergence.However, its convergence analysis is limited to the triangular meshes.The virtual finite element method [25] was also developed for solving two-dimensional Maxwell interface problems, in which each interface element is divided into subtriangles and the local space on the interface element is defined by piecewise Nédélec elements.However, its convergence analysis is also limited to the triangular meshes.The immersed WG method on triangular meshes was proposed in [59], and extended to polygonal meshes in [62].Compared to the lowest-order unfitted HHO method, the immersed WG method requires less restrictive mesh assumptions: the unfitted HHO method requires that two subregions of each interface element divided by the interface must contain a ball with radius comparable to the diameter of the element, while the immersed WG does not require such conditions.However, the immersed WG requires an additional regularity assumption: The Darcy velocity must be  1 on the entire domain.
In this paper, we define and analyze the lowest-order nonconforming immersed VEM for elliptic interface problems on unfitted polygonal meshes.Motivated by the conforming immersed VEM [26] and the nonconforming VEM [5], we define the virtual elements on each interface mesh element by the solutions of local interface problems with Neumann boundary conditions, and the elliptic projection is modified so that its range is the space of broken linear polynomials satisfying interface conditions, which is also used in the linear immersed FEMs (see, e.g., [52,55]).We derive optimal error estimates in the broken  1 -norm and  2 -norm under the standard regularity assumption that the solution is a piecewise  2 -function.Moreover, as in the immersed WG [62], the mesh assumptions in our scheme allow small cut elements.In addition, since our scheme is also a nonconforming method, the Darcy velocity can be recovered efficiently by casting a mixed formulation into the nonconforming method (see, e.g., [3,4,47,49,52]).We also note that, since there is an equivalence relation between the nonconforming VEM and the HHO method [34,54], we can reformulate our scheme in the context of HHO methods.However, it will be different from the unfitted HHO methods in [22,23], since these methods use the Nitsche-type formulation, while our method does not.
The rest of the paper is organized as follows.In Section 2, we introduce the model problem and mesh assumptions.In Section 3, we explain the nonconforming immersed VEM.In Section 4, we prove some approximation properties of the discrete spaces.In Section 5, we prove optimal error estimates of our scheme in the broken  1 -norm and the  2 -norm.Finally, we report several numerical tests that confirm the theoretical results in Section 6.

Preliminaries
We follow the standard notation of Sobolev spaces (see, e.g., [17,30]).For an integer  ≥ 0 and a subset  of R or R 2 , we denote by P  () the space of all polynomials of degree at most  on .For a subset  of R or R 2 , the indicator function on  is denoted by   .For a bounded measurable subset  of R or R 2 and  ∈  1 (), we denote by ()  the average of  on .
We consider the elliptic interface problem: with the jump conditions on the interface where the coefficient  is positive and piecewise constant on Ω ± , that is,   := | Ω  is constant for  = +, −.
for some generic positive constant  Ω .

Mesh assumptions
Let  ℎ be a decomposition (mesh) of Ω into polygonal elements  with maximum diameter ℎ.Let ℰ ℎ be the set of all edges in  ℎ .Let ℰ ∘ ℎ and ℰ  ℎ denote the set of all interior and boundary edges in  ℎ , respectively.For each  ∈  ℎ , let ℎ  and || be the diameter and the area of , respectively.For each  ∈ ℰ ℎ , we denote by || the length of .
An element  ∈  ℎ is called an interface element if the interface Γ passes through the interior of ; otherwise  is called a non-interface element.We denote by   ℎ and   ℎ the collections of all interface and non-interface elements in  ℎ , respectively.Analogously, an edge  ∈ ℰ ℎ is called an interface edge if Γ passes through the interior of ; otherwise  is called a non-interface edge.The collection of all interface edges and non-interface edges in ℰ ℎ are denoted by ℰ  ℎ and ℰ  ℎ , respectively.For each  ∈  ℎ , let ℰ   and ℰ   be the set of all interface and non-interface edges of , and let ℰ  := ℰ   ∪ ℰ   .We assume that ℎ is sufficiently small, and  ℎ satisfies the following regularity assumptions [5,7,21].Assumption 2.2.There exists  > 0 independent of ℎ such that (i) the decomposition  ℎ consists of a finite number of nonoverlapping polygonal elements; (ii) every element  ∈  ℎ is star-shaped with respect to a ball   with center   and radius ℎ  ; (iii) for each element  ∈  ℎ , all the edges in  have length larger than ℎ  ; (iv) the interface Γ meets the edges of each interface element at no more than two points; (v) the interface Γ meets each edge in ℰ ℎ at most once, except possibly it passes through two vertices.
Remark that the assumptions (iv) and (v) are reasonable for sufficiently small ℎ.As mentioned earlier, these assumptions allow small cut elements and do not require additional mesh procedures.
For each  ∈  ℎ , let   be the exterior unit normal vector along .For each  ∈ ℰ ℎ , let   be a unit normal vector of  with orientation fixed once and for all.Let  ∈ ℰ ∘ ℎ and let  1 and  2 be the polygons in  ℎ having  as a common edge.For  : Ω → R satisfying | 1 ∈  1 ( 1 ) and | 2 ∈  1 ( 2 ), we define the jump of  on  by For  ∈ ℰ  ℎ , we define []  := |  .We use the notations ∇ ℎ and ∇ ℎ • when the gradient and divergence operators are taken elementwise for piecewise smooth functions on  ℎ .

Nonconforming immersed virtual element method
In this section, we present the nonconforming immersed virtual element method for the elliptic interface problem (2.3).

Broken linear polynomials
We consider the space of piecewise linear polynomials satisfying the interface conditions on each interface element.Let  ∈  ℎ be an interface element.The broken polynomial space ̂︀ P 1 () is defined by It is easy to see that dim ̂︀ P 1 () = 3 and the following piecewise polynomials form a basis of ̂︀ P 1 (): and  = (− 2 ,  1 ).Since ̂︀ P 1 () ⊂  1 (), the space ∇ ̂︀ P 1 () is well-defined, and ∇ 2 and ∇ 3 form a basis of ∇ ̂︀ P 1 ().We set ̂︀ P 1 () := P 1 () for non-interface element  ∈  ℎ .

Nonconforming immersed virtual elements
We first define the local space on each element.For each interface element  ∈  ℎ , let where ̃︀ P 0 () is defined by ̃︀ P 0 () := {  + +   − : ,  ∈ R} for  ∈ ℰ  ℎ , and ̃︀ P 0 () := P 0 () for  ∈ ℰ  ℎ .Proceeding as in Lemma 3.1 in [5], it is easy to verify that ̂︀ P 1 () ⊂  ℎ () and the DOFs of  ℎ () can be chosen as follows: For each non-interface element  ∈  ℎ , the local space is defined as in the standard nonconforming VEM [5]: and its DOFs can be chosen as (3.1c).
Remark 3.1.If the local space  ℎ () on the interface element  defined as in the local space in the noninterface element, then the inclusion ̂︀ P 1 () ⊂  ℎ () and the interface condition [ ℎ /] Γ  ℎ = 0 does not hold.It may lead difficulties in analysis and extending our scheme to some other applications such as the nonhomogeneous interface problems.
With the above preparations, we state the nonconforming immersed virtual element method as follows: Find where the loading term ⟨, •⟩ is given by ⟨,  ℎ ⟩ := (, Π ∇  ℎ ) 0,Ω .Note that the well-posedness of the discrete problem (3.4) follows from Lemma 3.2.
Remark 3.3.The treatment of nonhomogeneous interface conditions is possible using such techniques as in [1,27].However, its analysis involves more technical issues.It is left for a future investigation.

Approximation properties on interface elements
In this section, we present some approximation properties of the broken linear polynomials and the immersed virtual elements on the interface elements.Note that some estimates were given in [62], but only the special case of the piecewise straight interface was considered.In this paper, we consider the curved interface.From now on, for ,  ≥ 0, we write   or   if there exists a constant  depending only on ,  and Γ (but independent of the location of the interface intersected with the mesh elements), and   if both   and   hold.
Using this lemma, we obtain the following projection error estimate.
Remark 4.4.Although it is not used in the following sections, one can prove the following interpolation error estimates: For  ∈   ℎ ,  ∈ ̃︀  2 Γ (), and

Error analysis
In this section, we derive the error estimate in the discrete energy norm for the scheme (3.4).We first compute the consistency error.
The following lemma shows that the discrete bilinear form  ℎ (•, •) is continuous on  ℎ (Ω) with respect to the  1 -seminorm.
We now present the error analysis in the energy-norm and the broken  1 -norm.
where the hidden constant also depends on  Ω in (2.4).

Numerical tests
In this section, we present some numerical tests for our proposed method.
Here the meshes in M3 are generated from PolyMesher [63].Some examples of the meshes are given in Figure 3.We compute the  1 -seminorm and  2 -norm errors given by where  ℎ is the solution of our scheme (3.4) (NC-IVEM).For all the examples, we plot the errors versus ℎ in Figures 4 to 6.We observe that the errors converge with optimal order, which is consistent with the theoretical result (see Thm. 5.4, Rem.5.5 and Thm.5.6).In particular, for Example 6.1, we also plot the results of the immersed WG method (IWG) [62] and the nonconforming immersed FEM (NC-IFEM) [46,52] on the uniform triangular meshes in Figure 7.In contrast with the NC-IVEM, both exhibit suboptimal convergence orders.Next, we compute the condition number for each example, where the meshes are fixed as M1.We plot the results in Figure 8.We observe that the condition number behaves as usual order (ℎ −2 ) for all the examples.
Here, the results indicated with the word "mod" will be explained later.We observe that the condition number seems to be proportional to  −1 .Further theoretical investigation is needed to explain such a phenomenon, which will be a subject of future work.In contrast, the errors remain bounded, which is consistent with our theoretical result that the error bounds are independent of the location of the interface intersected with the mesh elements (see Thm. 5.4, Rem.5.5 and Thm.5.6).
As observed above, the stiffness matrix becomes ill-conditioned in the presence of small-cut edges.However, we can improve the condition number by a simple modification of the discrete interface Γ ℎ like committing a variational crime in FEM: If  is an interface edge with | − | or | + | ℎ 2 , then we relocate the intersection point of  and Γ ℎ to the closest vertex of , and regard  as a non-interface edge (see Fig. 12).This procedure allows us to avoid extremely small-cut edges.The errors of the modified scheme unchanged (see Fig. 11 with "mod").
We next compute the condition number of the modified scheme.The results are reported in Figure 10 and are indicated with "mod."We observe that the condition number is uniformly bounded with respect to .

Test case 3: effect of small-cut cells
We also investigate the effect of small-cut cells, as in [35].Consider the problem with Ω = (0, 1)  13).In Figures 14 and 15 we plot the condition number and the errors versus  for each pair ( + ,  − ).Similar with the previous numerical test, we observe that the condition number seems to be proportional to  −1 , and the errors remain bounded.

Conclusion
We present the lowest-order nonconforming immersed VEM for elliptic interface problems with unfitted polygonal meshes.The local shape functions on the interface elements are defined as solutions of local interface problems with Neumann boundary conditions.We prove that our scheme achieves an optimal convergence rates in the broken  1 -norm and  2 -norm, under the piecewise  2 -regularity assumption.Some numerical tests are carried out to verify the theoretical results.
and let  ℎ Γ and  ℎ Γ be the unit normal and tangential vectors along Γ ℎ .

Figure 9 .
Figure 9. Test case 2: the interface (red line) and the mesh (black lines).

Figure 10 .
Figure 10.Test case 2: the comparison of condition numbers versus  with and without the discrete interface modification procedure.

Figure 11 .
Figure 11.Test case 2: the comparison of errors versus  with and without the discrete interface modification procedure.

Figure 13 .
Figure 13.Test case 3: the interface (red line) and the mesh (black lines).

Figure 14 .
Figure 14.Test case 3: the comparison of condition numbers versus .

Figure 15 .
Figure 15.Test case 3: the comparison of errors versus .