Maxwell interface problems widely appear in a large variety of science and engineering applications. In this article, we propose a virtual element method (VEM) to solve a two-dimensional (2D) H(curl; ⌦)-elliptic interface problem that is originated from Maxwell equations. One distinct advantage of the proposed method is its flexibility on the mesh generation to cater the interface. The mesh for computation is obtained from a background unfitted mesh by cutting interface elements into triangles and quadrilaterals, and the optimal convergence order is guaranteed independent of the mesh anisotropy.
To describe the idea, we let ⌦ ✓ Rfoot_0 be a bounded domain and let ✓ ⌦ be a smooth interface curve, as illustrated by the left plot in Figure 2.1. The interface cuts ⌦ into two subdomains ⌦ ± occupied by media with di↵erent magnetic and electric properties. We consider the following H(curl; ⌦)-elliptic interface problem for the electric field u :
with f 2 L 2 (⌦), subject to the Dirichlet boundary condition:
where the operator curl is for vector functions v = (v 1 , v 2 ) | such that curl u = @ x1 v 2 @ x2 v 1 while curl is for scalar functions v such that curl v = (@ x2 v, @ x1 v)
| with " | " denoting the transpose herein. The coe cients ↵ = ↵ ± and = ± in ⌦ ± are assumed to be positive piecewise constant functions of which the locations of the discontinuity align with one another. Moreover, we consider the following jump conditions at the interface :
[↵curl u] := ↵ curl u + ↵ + curl u = 0, (
where t denotes a tangential vector to . The condition (1.1c) is due to the tangential continuity of H(curl) functions and (1.1d) is from the fact that H(curl) is isomorphic to H 1 in 2D. The interface model (1.1) arises from each time step in a stable time-marching scheme for the eddy current computation of Maxwell equations [2,4]. In this model, ↵ denotes the magnetic permeability and ⇠ /4t is the scaling of the conductivity by the time-marching step size 4t.
For Maxwell equations, H(curl)-conforming Nédélec finite element spaces are widely used [17,25,37]. As for interface problems, the authors in [26] analyze the standard finite element methods (FEMs) for H(curl)-elliptic equations. The semidiscrete analysis for Maxwell interface problems with low regularity is provided in [43]. In addition, due to the potentially low regularity, there are many works focusing on developing a posteriori error estimators and adaptive FEM [8,16,22].
Numerical methods for solving interface problems based on unfitted meshes are attractive since they circumvent the burden of generating high-quality interfacefitted meshes which may be time-consuming for three dimensions (3D) or for moving interface problems. There have been a lot of works in this field on solving H 1 -elliptic interface problems, see [7,23,33] and the reference therein. However, there are much fewer works on solving Maxwell interface problems. Typical examples include matched interface and boundary (MIB) formulation [44], adaptive FEMs [15], and non-matching mesh methods [10,11,14]. Recently, a penalty method is developed in [35] requiring higher regularity.
A main di culty for solving H(curl) problems comes from low regularity of the exact solution, even for our assumption u 2 H 1 (curl; ⌦) := {u 2 H 1 (⌦), curl u 2 the polygonal mesh that satisfies a maximum angle condition [1]. Locally on each polygonal element, the new virtual functions can be also considered as discrete harmonic extensions according to the boundary conditions (degrees of freedom), while the usual virtual functions in [19,20] may be considered as continuous extensions. As the key advantage of using this new space, we are able to establish local Poincaré-type inequalities and optimal approximation capabilities for a large class of polygonal-shape elements all independent of element anisotropy which are the crucial intermediate results toward the final optimal error bound. A related work is [16] that constructs sub-meshes on interface elements of a background unfitted mesh for computation. One essential di↵erence between the proposed method and [16] is that the virtual mesh and space are only used for analysis in our work, while the computation procedure is the same as the usual VEM with the lowest order. The convergence is guaranteed independent of the element shape provided that the background mesh is shape-regular.
This article consists of 5 additional sections. In the next section, we introduce the background unfitted mesh and the fitted mesh according to interface geometry. In Section 3, we describe virtual spaces and projection operators. In Section 4, we present the numerical scheme and derive the error equation. In Section 5, we estimate the interpolation errors. In Section 6, we show that the convergence is of optimal order. In the last section, some numerical examples are presented to verify the theoretical estimates.
In this section, we first describe an unfitted background triangular mesh and then locally partition it into a fitted mesh used for computation. We then introduce Sobolev spaces for H(curl)-interface problems. Although the triangular and quadrilateral shape is the focus of this work, we highlight that most key results are actually established and presented for more general polygonal element shapes.
Consider an interface-independent shape-regular triangular mesh of the domain ⌦. Note that this mesh could be simply taken as a highly-structured mesh due to the interface independence. We shall call it a background mesh, and denoted it by T B h . Another example of T B h is a uniform Cartesian grid which is widely used for unfitted mesh methods. The proposed analysis approach can be easily adapt to this grid. If a triangular element in T B h intersects the interface, then it is called an interface element. The collection of interface elements is denoted as T Bi h . The remaining elements are called non-interface elements. For this background mesh, we further make the following assumptions: connecting the two intersection points. Then all these connected small segments, defined as h , form a piecewise linear approximation of the true interface . In addition, by the assumption (A), each triangular interface element K is cut by K h into a triangular and a quadrilateral subelement from which an interface fitted mesh can be generated. We denote this fitted mesh by T h . The collections of quadrilateral and triangular elements in T h resulted by the interface-cutting are denoted by T q h and T t h , respectively. Those elements all have one edge aligning with the interface approximately, and thus are also called interface elements in T h . Clearly, there holds
h are only di↵erent on the interface elements. Furthermore, an interface element K is assumed to be cut into K h and K + h by K h , and the mismatch portion, with K ± := K \ ⌦ ± cut from the original interface, is denoted by K int , indicated by the shaded region in the left plot of Figure 2.2.
In addition, for each interface edge K h , we assume there is a shape regular triangle B K h ✓ ⌦ with the base K h and a height O(h K ). Here B K h is not required to align with the elements in the mesh. Further assume all the B K h have finite overlapping. Note that for the considered background regular triangular mesh, for each interface element K, this B K h certainly exists and can be further shown to be contained in K.
Moreover we assume the interface is well-resolved by the mesh, and it can be quantitatively described in terms of the following lemma [23].
Lemma 2.1. Suppose the mesh is su ciently fine such that h < h 0 for some valve value h 0 , on each interface element K 2 T Bi h , there exist a constant C independent of the interface location inside K and h K such that for every point x 2 \ K with its orthogonal projection x ? onto K h , dist(x, x ? ) Ch 2 K .
(2.1) Fig. 2.2: Left: an interface element K 2 T B h is cut to a triangular element K t and a quadrilateral element K q of which both are in T h . Right: the quadrilateral element is further cut into two triangular elements K q 1 and K q 2 .
The explicit dependence of h 0 on the curvature of the interface can be found in [23]. An adaptively-generated background mesh can be found in [42] capturing large curvature of interface curve. Now following [26,34], we introduce the -strip:
By the estimate (2.1), we have
with the constant C only depending on the interface. Furthermore, from [26,34] we can control the L 2 -norm in the -strip by the width of the strip and thus obtain first order convergence when = O(h 2 ).
Lemma 2.2. It holds for any
We next introduce some major Sobolev spaces used throughout this article. For each subdomain ! ✓ ⌦, we let H s (!) and H s (!), s 0, be the standard scalar and 2D vector Hilbert spaces on ! where specifically H 0 (!) = L 2 (!) and H 0 (!) = L 2 (!). In addition, for s 0, we let
(2.5)
Similarly, we introduce H s (div; !) as the counterpart of H s (curl; !) with the divergence operator. If ! \ 6 = ;, then ! ± = ! \ ⌦ ± , and H s (curl; ! [ ! + ) denotes the space of functions piecewisely defined in H s (curl; ! ± ). For these spaces, we can define their subspaces H s 0 (!), H s 0 (!), and H s 0 (curl; !) with the zero trace on @!. Also let (•, •) ! be the standard L 2 inner product on !.
When the interface is smooth, the solution to the problem is expected to have an H 1 (curl; ⌦ ± ) regularity ( [18,30]). The fundamental H 1 (curl; ⌦)-extension operator established by Hiptmair, Li and Zou in [26] (Theorem 3.4 and Corollary 3.5) will be used in analysis.
Theorem 2.1 (Theorem 3.4 and Corollary 3.5 in [26]). There exist two bounded linear operators
such that for each u 2 H 1 (curl; ⌦ ± ):
1.
with the constant C E only depending on ⌦ and .
Using these two special extension operators, we can define u ± E = E ± curl u ± which are the keys in the analysis later. Finally, throughout this article, for simplicity we shall use . to denote a • • • C • • • with a generic constant independent of mesh size and interface location relative to the mesh.
In this section, we shall introduce a virtual element space using the lowest order Nédélec element on a virtual triangulation which is obtained by refinement of the background mesh.
In the proposed method, triangular elements and quadrilateral elements in T h are treated di↵erently. To avoid confusion, in this section we shall usually use K t and K q to denote triangular and quadrilateral elements in T t h and T q h , respectively, while K denotes an interface element in the background mesh T B h or general elements in T h if there is no need to distinguish their shape. In addition, for simplicity's sake, we always use h K as the diameter of elements K, K q and K t .
For any element or edge ! in T h , we let P k (!) be the k-th degree polynomial space defined on !. Given a triangle K, we shall consider the first family Nédélec element of the lowest degree [39] as the underling approximation space:
Then (3.1) will be used on all the non-interface elements of T B h (or T h ) as well as the triangular interface elements K t . But for K q we need to employ a virtual element space. Let us first discuss the definition on a general polygon P : Ṽh (P ) = {v h 2 H(curl; P ) \ H(div; P ) : v h • t e 2 P 0 (e), 8e ⇢ @P,
This is exactly the one introduced in [20,19] with the lowest degree and N D h (P ) ⇢ Ṽh (P ). Similar to the ones for (3.1) defined on triangles, the local degrees of freedom (d.o.f.) for (3.2) are v h e • t e , e ⇢ @P.
(3.
3)
It has been shown in [20,19] that the functions in (3.2) can be uniquely determined by d.o.f. (3.3). For the present situation, Ṽh (K q ) (P = K q ) is the space used for discretization on K q . As the shape of elements K q could be very anisotropic, a robust norm equivalence and interpolation error estimate is hard to establish in Ṽh (K q ). To address this issue, we shall introduce an auxiliary triangulation of ⌦ (a virtual mesh), and construct an auxiliary H(curl)-conforming space associated with this mesh. Given an interface element K 2 T Bi h , with an quadrilateral subelement K q 2 T q h , then the local auxiliary mesh is formed by a Delaunay triangulation of K q : connecting the diagonal s.t. the sum of angles opposing to the diagonal is less than or equal to ⇡; see the right plot in Figure 2.2 for an illustration. Although it may contain anisotropic triangles, each triangle from this new partition satisfies the maximum angle condition (Lemma 3.1), which is key to robust interpolation error estimates (see Section 3.2). A similar result is proven in [13] for Cartesian grids. Lemma 3.1. Let K be a shape regular triangle, i.e., there exist 0 < ✓ min ✓ max < ⇡ such that every angle ✓ in K satisfies ✓ min ✓ ✓ max , then every triangle in the auxiliary Delaunay triangulation on K described above satisfies the maximum angle condition, i.e, every ✓ in the auxiliary triangulation is bounded above by ✓max max{⇡ ✓ min , ✓ max }.
Proof. Without loss of generality, we consider the triangle in the right plot of Figure 2.2 for illustration where the left and right cutting points are D and E, respectively. We shall bound angles of three triangles:
If the angle is one of (or part of) the angles of A 1 A 2 A 3 , then it is bounded by ✓ max . If the triangle contains one of the angles of A 1 A 2 A 3 , as the sum of three angles is ⇡, we conclude other angles are bounded by ⇡ ✓ min . So we only focus on angles of the triangle DA 2 E. Now use the Delaunay property,
The discussion for the special case on a quadrilateral K q is postponed to Section 5, and the results on a general polygon P are the focus in the rest of this section. To be able to establish a robust analysis of the approximation capabilities of V h (P ) below and the stability of the discretization, P is assumed to admit a triangulation T h (P ) with no additional interior vertices added, i.e., the collection of edges E h (P ) in T h (P ) are solely formed by the vertices on @P , and this triangulation satisfies (P1) the maximum angle condition for each triangle in T h (P ); (P2) no-short-interior-edge condition: h P . |e| for every interior edge e;
(P3) star-convexity: there exists x := (x 1 , x2 ) 2 P such that xy ✓ P , 8y 2 @P .
We then define an auxiliary VEM space
Indeed we can show the VEM space defined by (3.4) shares the same d.o.f. of (3.3).
Lemma 3.2. Let P be a simple polygon, then the d.o.f. v h •t e , e ✓ @P are unisolvent on the space V h (P ).
Proof. Let us consider the following well-posed problem: for any given boundary conditions v h • t e , e 2 @P , find (v h , h ) satisfying
where S h,0 (T h (P )) is the piecewise linear Lagrange finite element space. Since there is no internal vertex, S h,0 (T h (P )) is a trivial space, and (3.5) reduces to
The fact that curl v h is a piecewise constant on T h (P ) and integration by parts show X e2E h (P )
Therefore, curl v h must be a single constant on all elements in T h (P ). Namely, the solution space of (3.5) is V h (P ), and the unisolvence follows from the homogeneous boundary condition yielding the zero solution.
Remark 3.1. The proof of Lemma 3.2 basically shows that v h 2 V h (P ) satisfies div h v h = 0 together with (3.6), where div h is the element-wise div operator. Namely, the functions in the VEM space (3.8) can be treated as a discrete harmonic extension of the boundary conditions v h • t on @P while the original virtual space (3.2) is a continuous extension with a pointwise constraint div v h = 0. Therefore functions in Ṽh (P ) may not be polynomials while V h (P ) has piecewise polynomial vectors for which the error estimates is relatively easy to establish.
By Lemma 3.1, it is clearly that K q = K q 1 [ K q 2 induces such a triangulation satisfying (P1)-(P3). Meanwhile, we would like to remark that the aforementioned setting and the forthcoming analysis in this paper can be easily extended to the case when T B h is a uniform Cartesian grid, on which the interface elements consist either trapezoids, or triangle-pentagon satisfying (P1)-(P3). In the subsequent analysis involving K q , the space
All these triangles on interface elements form a triangulation resolving the interface, and the global H(curl)-conforming space is defined as
(3.9)
As we assume the background mesh T B h is shape regular, the maximum angle condition holds uniformly for the auxiliary mesh T h
For a general polygon P and V h (P ), the constant curl v h in P can be computed by d.o.f. as
With curl v h and v h • t known, its L 2 -projection can be computed following [20,19]. On any elements or edges ! ✓ ⌦ we define the local
which is indeed computable according to the d.o.f. of (3.2) [19, Remark 3]. For readers' sake, we recall the procedure here: for each p = (p 1 , p 2 ) | 2 [P 0 (P )] 2 , there exists h = p 2 (x 1 x1 ) + p 1 (x 2 x2 ) 2 P 1 (P ), such that curl h = p, where (x 1 , x2 ) is the point in the star-convexity assumption (P3), Therefore
(3.12)
As v h • t is given as d.o.f., and curl v h is constant, we get
in which the integration on @P is with respect to ds(x 1 , x 2 ). Due to the d.o.f. imposed on edges, we can define the interpolation
We note that if P is a triangle, I P reduces exactly to the usual edge interpolation operator, and the special one is for other general polygons such as quadrilateral elements K q where shape functions are from the virtual space V h (P ) in (3.8). Using integration by parts, we get
Namely curl I P u is the L 2 -projection of curl u to the space of constants. Moreover the interpolation I P : Ṽh (P ) ! V h (P ) serves as a bijective mapping which also preserves curl values and the L 2 projection onto [P 0 (P )] 2 as both Ṽh (P ) and V h (P ) share the same d.o.f. . For the considered mesh T h , taking P = K 2 T h , we have Ṽh and V h lead to the same numerical scheme but the analysis based on
V h can exploit more existing tools built for simplicial finite elements. Finally, a global interpolant u I is formed by gluing these local interpolations together, for which certain modification must be introduced on the interface edges forming h (see Section 5.2).
In the rest of this section, we present some estimates which show the convenience in analysis of opting for the space V h (P ). For a triangle with vertices a i , let ✓ i be the angle at vertex a i and e i be the edge opposite to a i , for i = 1, 2, 3.
Lemma 3.3. The following identity holds for any linear h on a triangle T :
where R T is the circumradius of T and t i is a unit tangential vector of e i .
Proof. Denote by i := h (a i ) for i = 1, 2, 3. The cotangent formula [12] reads
Then the law of sines and
We now prove the following Poincaré-type inequality which is the key for the analysis on anisotropic meshes. Lemma 3.4. Let P be a simple polygon satisfying (P1)-(P3), then
(3.16)
Proof. Define an auxiliary function
, where (x 1 , x2 ) is the point in the star-convexity condition in (P3). It is clearly that
In addition, for every edge e 2 @P , ( (
• t| e yields the height l e of e in the triangle formed by e and (x 1 , x2 ), thus we have
Together with the star-convexity condition, we have
Note that curl(w h v h ) = 0, then by a standard argument of the conforming exact sequence, there exists a continuous piecewise linear finite element function h s.t.
v h w h = r h . Applying Lemma 3.3, together with the maximum angle condition in (P1), we get the estimate
We then control the norm contribution from an interior edge e. Since P is simply connected, any interior edge e divides P into two parts. Choose the part with less boundary edges and denote it by P e . Note that R @Pe r h • t ds = 0, consequently by r h • t being a constant on each edge on @P e , we have an identity decomposing @P e = (@P e \ @P ) [ e,
and thus
Then from (3.19) and the condition (P2) we can get
⌘ , with constant depending on the number of vertices of P but not the shape regularity of P . Finally, the desired estimate (3.16) follows from the triangle inequality and estimates (3.17)-(3.18).
In this section, we describe the proposed virtual element formulation and derive an error bound. We start with the standard weak formulation: find u 2 H 0 (curl, ⌦) such that
We emphasize that the local "virtual" element space (3.8) and the global one (3.9) is right away a computable space, readily used for the discretization, unlike (3.2). The d.o.f. on the diagonal edge can be determined by solving (3.6) explicitly, and a set of modified harmonic bases on boundary edges can be obtained and used in computation. As a result, the standard Galerkin formulation is computable without referring to the VEM framework of a projection-stabilization split:
where ↵ h and h are the modification of ↵ and according to the linearly approximated interface h . No projection operator is required since all the shape functions are computable. However, this approach will introduce an extra partition which becomes inecient especially in 3D. Instead, we shall treat henceforth the interface part of T h as a virtual mesh only appearing in analysis not computation, whereas this associates the meaning of "virtual" in V h . Its approximation capabilities will be discussed in Section 5.1 based on the maximum angle condition.
Using the L 2 -projection (3.11), we define a bilinear form
where the operator ⇧ h is taken as ⇧ K q if K = K q 2 T q h , and the identity operator otherwise. The stabilization S K (u h , v h ) is defined element-wisely only on K q 2 T q h , i.e., the quadrilateral subelements of the interface elements K 2 T Bi h :
with a parameter K independent of the mesh size and specified later. Note that the motivation of this stabilization term comes from the approximation of ( h (u h ⇧ h u h ), v h ⇧ h v h ) K q and thus suggests the scaling h K in (4.4).
At last, the proposed VEM discretization is to find u h 2 V h such that
As mentioned in Section 3, some elements could be extremely anisotropic, and the commonly used norm equivalence in the VEM framework may not be applicable. Following the approach in [9], we shall work on an induced norm on V h by the bilinear form in (4.3) (Lemma 4.1) which is weaker than the original graph norm:
Proof. Suppose |||v h ||| h = 0, then clearly v h = 0 on all triangular elements in T h . So we only need to consider K q 2 T q h . Indeed, ⇧ K q v ⌘ 0 on K q and (v h ⇧ K q v h )•t vanishing on @K q implies v h •t = 0 on @K q . Due to the unisolvence, we have v h = 0 on K q which finishes the proof.
In the following main theorem, we derive an error equation to (4.5) to demonstrate how the VEM framework can in a novel manner overcome the di culties of the non-conformity issue aforementioned in the introduction for other DG-based approaches. To reinstate the optimal rate of convergence, we need to further assume that the source term f bears certain extra local regularity. First, the error is decomposed and the error equation is for ⌘ h :
where u I 2 V h is an arbitrary function in VEM space.
) is the solution to (4.1) and let u I 2 V h be an arbitrary function in VEM space, then for
Proof. We have
. (4.9)
For (I), on all the triangular elements in T h , ⇧ h reduce to identity operators, so ⇧ h ⌘ h ⌘ h simply vanishes. On a quadrilateral element K q 2 T q h , by the definition under (4.3),
For the term (II) in (4.9), using ↵ curl curl u + u = f , we have
For (IIa), since ⌘ h is in the conforming auxiliary space V h in (3.9), using the integration by parts, the continuity condition of the original PDE, and the curl condition in (3.8) we immediately have
(4.12)
Instructions for Typing Manuscripts (Paper's Title) 15
For (IIb), on a triangular element K t , we note that
On a quadrilateral element K q , by (3.11) we have
Combining (4.13) and (4.14), we have
In addition, for the stabilization term (IIc), by (⌘
• t e 2 P 0 (e) on each e ✓ @K q and the definition of the interpolant in (3.14), we have
Finally, putting the estimates in (4.10)-(4.16) to (4.9) yields the following bound
⌘ .
(4.17)
To bound k⇧ K q ⌘ h ⌘ h k L 2 (K q ) on quadrilateral elements, using Lemma 3.4 yields the following estimate
. Putting the estimate above into (4.17) and canceling one |||⌘ h ||| on both sides yield the desired result.
In this section, we estimate the interpolation errors and projection errors of virtual element spaces. Given any triangle T , the interpolation in (3.14) exactly becomes the canonical edge interpolation [38]. If T is further assumed to be shape regular, then the following standard optimal approximation capability holds:
(5.1)
Due to the assumption of the interface being smooth, we note that certain elements in T t h may inevitably have high aspect ratio in the process of mesh refining, which results that the commonly assumed shape regularity does not hold anymore. Consequently, the standard results about the approximation results of the edge interpolation (5.1) cannot be directly applied. However, since maximum angles of triangles in the auxiliary triangulation around the interface are uniformly bounded if the background mesh is shape regular, the interpolation error estimates can nevertheless be established based on the maximum angle condition. The interpolation estimates based on the maximum angle condition have been long studied for Lagrange elements [3,32], Raviart-Thomas elements [1,6,36], and 3D Nédélec elements [6].
Lemma 5.1 (the same argument as in [1]). Given any triangle T , let ✓ T be the maximum angle of T , then
2)
The results above can be directly applied to estimate the interpolation errors of the virtual space V h (K q ) on K q 2 T q h . Again we present in a more general setting.
Lemma 5.2. Let P be a simple polygon satisfying (P1)-(P2) and let I P u be the edge interpolation to V h (P ) defined in (3.14). Then ku I P uk H(curl;P ) . h P kuk H 1 (curl;Conv(P )) , u 2 H 1 (curl; Conv(P )). ( 5.3)
Proof. The estimate for the semi-curl norm is standard since curl I P u is the L 2 projection of curl u on P . Then
where the last step is the Poincaré inequality over convex domains [40]. Let I h be the edge interpolation to V h (T h (P )), i.e., the standard edge finite element space on mesh T h (P ). By the maximum angle condition in (P1) and Lemma 5.1, we have ku I h uk L 2 (P ) . h P kuk H 1 (curl;P ) . Then it su ces to estimate the di↵erence kI P u I h uk L 2 (K) on each triangle K 2 T h (P ). We apply Lemma 3.4 on each K to get
As (I P u I h u) • t = 0 for e ⇢ @P , we only consider an interior edge e. Since P is simple, any interior edge e divides P into two parts. Choose the part with less
Using the triangle inequality, together with the estimates for curl(u I P u) and curl(u I h u) , we conclude for any K 2 T h (P )
(5.4)
The desired result (5.3) then follows from the triangle inequality.
In the interpolation error estimate, locally a norm k curl uk H 1 (K) will be used. When K is an interface element, in general curl u 6 2 H 1 (K) but in H 1 (K + [ K ). Instead we will use the fact curl u ± E 2 H 1 (K) and define the interpolation by the tangential components of either u + E or u E , where which extension to use depends on the measure of K and K + . Note that, in the present situation, since both the triangular elements in T t h and the quadrilateral elements in T q h may have high aspect ratio, the modification in [26] may not be suitable on anisotropic meshes with interface being present. Therefore, we shall employ a di↵erent interface-aware interpolation.
In the following discussion, we only present the results for the elements in the mesh T h due to the technical treatment for the interface. But we emphasize that most of the results can be generalized based on the estimate of the interpolation errors on general polygons above. We shall use K to denote an interface element in
h that is cut into K h and K + h by the edge K h , and without loss of generality, we assume K h 2 T t h and K + h 2 T q h . Recall that K int is the portion sandwiched between and K h , and we further define
namely the mismatching subregions of K ± h as shown in Figure 5.1. Let E K be the collection of edges of K h and K + h but excluding the edge K h . We define a modified interpolation operator ĨK on
By such a definition, we can always keep the interpolation as the standard one on the subelement with smaller size. So when estimation on the mismatch portion is needed, such as (5.15), the element size appearing on the denominator will be always larger than |K|/2 such that the overall estimate can be controlled. This consideration serves as our key motivation to make this modification. For simplicity, we denote u I by the global interpolant such that
In addition, we use I
E to denote the canonical interpolation on K ± h for Sobolev extensions u ± E . We emphasize that the modified ĨK serves the purpose for the error analysis and is not needed in actual computation.
The following two lemmas are presented for general polygons. So we temporarily let K be an interface polygon, and the notation K , K h and K int are all defined in the same manner as their counterparts for triangular interface elements. For the subelement with larger size, inevitably there is a mismatch on K h , so these results are essential. In the following discussion, with slightly abuse of the notations, we denote hK = | K h | which might be much smaller than h K (see Fig. 5.1 (left)), and
with trivial generalization to other Sobolev norms. Fig. 5.1: Left: the triangular interface is the smaller one. Right: the quadrilateral element is the smaller one.
h is assumed to be part in ⌦ + and part in ⌦ . As a result, the line integral R K h u•t ds has part of the integrand being u + E •t while the other being u E •t. Their di↵erence appears one of the key terms to bound the error of the modified interpolant (5.5), as one adheres to one extension in defining the interpolation.
Instructions for Typing Manuscripts (Paper's Title) 19
Proof. Applying integration by parts on K int and using the jump condition in (1.1c), we obtain
which yields (5.6) since |K int | . hK h 2 K by (2.1). The result above can be used to derive the following trace inequality. Recall that there exists a shape regular triangle B K h ✓ ⌦ with the base K h and a height O(h K ) by Assumption (B) for all interface elements. Lemma 5.4. Let u 2 H 1 (curl; ⌦ [ ⌦ + ). Given an interface polygon K with K h , there holds
Proof. Apply the L 2 -projection on K h to obtain k(u
.
(
Since t is a constant vector and K h with B K h satisfies the height condition, by the trace inequality [9, Lemma 6.3] and Poincaré inequality with average zero on a boundary edge [9,Lemma 6.11], we have
. h
(5.9)
For (II), by Lemma 5.3, we have
(5.10) Putting (5.9) and (5.10) back into (5.8) finishes the proof.
Now we proceed to estimate the interpolation errors u u I on interface elements for the modified interpolation. Lemma 5.5. Let u 2 H 1 (curl; ⌦ [ ⌦ + ). Given each interface element K 2 T Bi h , there holds
(5.11)
Without loss of generality, we focus the proof on K h as the estimate on the other part follows the result on K h using a similar argument as the one in Lemma 5.2. By the triangle inequality, we have
The first term (I) can be bounded by
and the second term (II) directly follows from Lemma 5.
. 1
To estimate the L 2 -norm, we use inequality (3.16) in Lemma 3.4 to conclude
(5.15)
Lastly, using Lemma 5.4 and the bound of kcurl w h k L 2 (K h ) finishes the proof.
The estimates on non-interface elements in the background mesh T B h are standard. These estimates together with the Sobolev inequality in Lemma 2.2 and Theorem 2.1 on the extension yield the global interpolation estimate.
Theorem 5.1. Let u 2 H 1 (curl; ⌦ [ ⌦ + ), then there holds
(5.16)
Proof. For non-interface elements, the estimate is standard as well. For interface element K, we then use Lemma 5.5:
, in which the second step we use the fact h is uniform Lipschitz so that the overlapping portions of triangles B K h for every interface element K are bounded. The desired estimate follows from Theorem 2.1 and estimate (2.4).
In this subsection, we move back to the mesh T h consisting of triangular and quadrilateral elements cut from the background triangular mesh. On the quadrilateral K + h , a stabilization term is present. In this section such terms appearing in the error bound are estimated, including
The main di culty is on the second term above. Note that the common and natural approach to estimate the edge terms is to apply the trace inequality, which indeed works for the edges A 1 D and A 2 E due to the corresponding O(h K ) height within the triangle. However, the major di culty arises for edges like A 1 A 2 and K h = DE, due to a possibly degenerating height. The core idea of our approach is to employ a constructive proof, without relying on the trace inequality, to control the edge terms by using ⇧ K q (u u I ) • t being a constant for lifting and applying the definition of projection (3.12). In the coming proofs, ⇠ h := u u I for simplicity.
Proof. First, we have on each edge e 6 = K h ✓ @K + h , R e ⇠ h ds = 0, and thus
where we have used the fact that e is one part of an edge of the regular element K such that the trace inequality can be applied on this edge and K. On K h , by the triangle inequality, we have
For the first term in (5.19), note that
Proof. For simplicity, we assume A 1 is at the origin, K is contained in the first quadrant, and the edge A 1 A 2 aligns with the x-axis having a tangential vector (1, 0) | . Let e be an edge of K + h with the unit tangential vector t e . If e = A 1 D or A 2 E, since the height within K + h with respect to these two edges cannot degenerate, a simple scaling directly leads to k⇧
.
(
of which the estimate follows from Lemma 5.5. For (II), on each e 0 ✓ @K + h using (5.22) Lastly, we arrive at
(5.27)
Proof. Let us decompose the error into
(5.28)
Here the estimate of the second term is similar to the one in Lemma 5.7. Therefore, we only need to estimate the first term in (5.28) which is further decomposed into
.
(5.29)
Note that (I) is only non-zero on K h of which the estimate follows form Lemma 5.4. For (II), if e ✓ @K + h is A 1 D or A 2 E, i.e., it has an O(h K ) height within K + h . Then we apply the trace inequality [9, Lemma 6.3] and the approximation result of the L 2 projection to obtain (II) . h
. If e ✓ @K + h is A 1 A 2 or DE where the corresponding height may become degenerate, we first apply the trace inequality on the whole shape-regular element K, and then apply the Poincaré inequality [9,Lemma 5.3], to obtain (II) . h
. Combining the estimates above finishes the proof.
In this section, based on the previous results, we estimate the convergence order of the solution errors. In particular, we need to estimate each term in the error bound (6.2). Our main task is to estimate those terms on quadrilateral elements. In the following discussion, we still keep our notation that K + h 2 T q h will be the quadrilateral subelement associated with each interface element K 2 T Bi h . Theorem 6.1 (An a priori convergence result for VEM). Under the same assumption of Theorem 4.1, let u 2 H 1 (curl; ⌦ [ ⌦ + ) and let the background mesh T B h satisfy the assumptions (A) and (B), then the solution u h to the VEM scheme (4.5) admits the error estimates ku u h k H(curl;⌦) . hkuk H 1 (curl;⌦ + [⌦ ) + h
Recall that in Theorem 4.1, we have obtained
Then the estimate follows from Lemma 5.8 and Theorem 5.1, and applying a simple triangle inequality to the last term.
In this section, we present a group of numerical experiments to validate the previous estimates. Let the computation domain be ⌦ = ( 1, 1) ⇥ ( 1, 1), and background mesh be generated by triangulating an N ⇥N Cartesian mesh by cutting each square into two triangles along its diagonal. We highlight that the proposed method can be used on any other regular background triangular meshes. A circular interface { : x 2 + y 2 = r 2 1 } cuts ⌦ into the inside subdomain ⌦ and the outside subdomain ⌦ + . We consider the example in [24,34] that the exact solution is
where the boundary conditions and the right hand side f are calculated accordingly. We employ the parameters k 2 = 20, k 1 = k 2 (r 2 2 r 2 1 ) with r 1 = ⇡/5 and r 2 = 1, and fix ↵ = = 1 with varying ↵ + = 10 or 100 and + = 10 or 100. For simplicity, we define the errors e 0 = ku u h k L 2 (⌦) and e 1 = k curl(u u h )k L 2 (⌦) .
(7.2)
Authors' Names
Authors' Names
Authors' Names