AN ADAPTIVE EDGE FINITE ELEMENT METHOD FOR THE MAXWELL’S EQUATIONS IN METAMATERIALS

. In this paper, we study an adaptive edge ﬁnite element method for time-harmonic Maxwell’s equations in metamaterials. A-posteriori error estimators based on the recovery type and residual type are proposed, respec-tively. Based on our a-posteriori error estimators, the adaptive edge ﬁnite element method is designed and applied to simulate the backward wave propagation, electromagnetic splitter, rotator, concentrator and cloak devices. Nu- merical examples are presented to illustrate the reliability and eﬃciency of the proposed a-posteriori error estimations for the adaptive method.

1. Introduction. In this paper, we consider the following two-dimensional time harmonic Maxwell's equations where Ω ⊂ R 2 is a bounded domain with Lipschitz boundary ∂Ω, E = [E 1 , E 2 ] t denotes electric field and F = [F 1 , F 2 ] t ∈ [L 2 (Ω)] 2 represents the given wave source term, where [L 2 (Ω)] 2 is denoted as the two-dimensional vector space of square integrable functions. The parameters µ and ε are the permeability and permittivity of the material medium, which may be functions of position. And k is the wavenumber.
is defined differently depending on the space dimension d. For d = 1, the curl of a scalar v is defined as the vector Solutions of many problems of interest require computation of phenomena occurring in a physical domain of infinite extent. In such an infinite domain, electromagnetic sources are all considered to be at the origin, and to scatter energy to infinity. A condition is required to prevent the nonphysical inverse process, where energy is transferred by radiation from infinity to the origin, and ensure uniqueness of the solution in an infinite domain. For the two-dimensional case, the following condition, known as the Silver-Müller radiation condition [5], must be satisfied lim r→∞ r 1/2 (∇ × E × n − ikE) = 0, r = |x|. ( The numerical approximations of the Maxwell's equations have been extensively investigated in the last many years, such as finite element method (FEM) [6,7,15,18,23,36], finite difference time domain method (FDTD) [14,30], interior penalty discontinuous Galerkin method (IPDG) [17], hybridizable discontinuous Galerkin method (HDG) [27], nodal discontinuous Galerkin method (NDG) [22] and multigrid method (MG) [11], etc. In recent years, metamaterials as a new type of material have attracted people's attention because of their potential applications in many fields [14,28], such as cloaking devices [2,12,24], sub-wavelength imaging and perfect lens. We mention in passing some related mathematical research on metamaterials in the literature [1,4,20,21], and some related finite element methods designed for simulating the electromagnetic wave in metamaterlals [33,34,35]. The advantage of adaptive finite element method in metamaterials, which will be considered in this paper, is to make possible, with reasonable computational cost, the local mesh adjustment of the computational domain according to a-posteriori error estimator.
The adaptive method can improve the efficiency of electromagnetic field calculation, and it is simple and and accurate [6,7,23]. In this paper, we develop the adaptive finite element method of Maxwell's equation in metamaterials. Because many electromagnetic problems are actually interface problems, we need to deal with them carefully. Some parts of the interface problems do not need to be refined near the interface, but they are still refined by using the usual indicators. So the more efficient a-posteriori error estimator should be considered without regard to the particularity of the examples and make the meshes match the images of the numerical solution better, which is the purpose of this article. In view of such problems, the recovery technique is used in each sub-area and the corresponding error is calculated. Errors of all sub-areas constitute errors of the domain of interest to guide the grid refinement.
The rest of this paper is organized as follows. In Section 2, we introduce the finite element scheme for the Maxwell's equations. A-posteriori error estimators are proposed, which include two types: one is the residual-based and the other is the recovery-based in Section 3. Afterwards in Section 4, several numerical examples are presented to compare the realization effect of the different posteriori error indicators. In Section 5, we give a conclusion.
2. Finite element scheme. Let L 2 (Ω) be the Hilbert space of square integrable functions on Ω equipped with the standard L 2 -norm: w = (w, w), where (·, ·) is standard L 2 -inner product. Consider the following Hilbert space where n is the unit outward normal vector to ∂Ω.
The Berenger PML [3] is applied to approximate the Silver-Müller radiation condition with the perfect electric conductor (PEC) boundary conditions The variational formulation of Eq.(1) reads: Given F ∈ [L 2 (Ω)] 2 , find E ∈ H 0 (curl; Ω) such that where Let T h = {K} be a regular triangle partition of the domain Ω with mesh size h, define the lowest order edge elements space [25] A basis for V h is given by the set of Nédélec's shape functions {φ i } ne i=1 with n e the number of triangle edges in the mesh, the simple formula for local shape functions is defined as , and associated with linear Lagrangian basis functions {λ i } 3 i=1 , where the explicit expressions for the basis λ i of the triangle K is given by with cyclic permutation of the indices {i, j, k} over {1, 2, 3} and |K| is the area of K.
It can be shown that edge finite elements guarantee the continuity of the tangential component across faces shared by adjacent triangles; they thus fit the continuity properties of the electric field. With the choice of discrete space V h , the curl-conforming finite element scheme of Eq.(4) takes the form as follows: 3. A-posteriori error estimations. Having computed a finite element solution, it is possible to obtain a-posteriori error estimations. These estimations accomplish two main goals. Firstly, they are able to quantify the actual error numerically. Secondly, they can be used to perform the adaptive mesh refinement, which solves the problem iteratively. The method uses a-posteriori error estimators to indicate where errors are particularly high, and more mesh intervals are then placed in those locations. This process Solve =⇒ Estimate =⇒ Mark =⇒ Refine can be repeated until a satisfactory error tolerance or result is reached. In this section, we propose a-posteriori error estimations including not only a residual type but also a recovery type. A Dörfler marking strategy [13] is used in the adaptive process.
3.1. The residual-based a-posteriori error estimation. The residual type aposteriori error estimator η r0 K l [15] for the element K l is defined as where E denotes the set of all interior edges, h K l = |K l | 1/2 is the the diameter of K l , and the local contributions R 1 K l , R 2 K l and J 1 E , J 2 E are defined, respectively, as

3.2.
The recovery-based a-posteriori error estimations. The recovery type a-posteriori error estimators η r1 K l [23] and η r2 K l are considered, respectively, as and for K l ∈ T h , here R denotes the recovery operator. The specific recovery algorithm see the Algorithm 1.

Algorithm 1
• For the interior edge E i , we select two elements K i1 , K i2 sharing it, i.e., where i = 1, 2. Take α 1 = 1, β 1 = 1 for Eq.(8) and α 2 = µ −1 , β 2 = ε for Eq.(9), respectively. • For the boundary edge E i , we find several interior edges that are closer to it and by linear interpolation, we can obtain R(α i ∇ × E) and R(β i E) at the midpoint x i of the edge E i .

Remark 1.
Recovery operator can be built in a variety of ways such as gradient recovery [19,26,37], functional recovery [29], flux recovery [8,9,10,31] and so on. Two quantities related to µ −1 ∇ × E and εE are recovered in the respective H(curl)-and H(Div)-conforming finite element spaces in [9,10]. It is different that we take average of flux to construct the flux recovery operator, which is easy to implement with superconvergence. (1) Assume that |µ −1 | < C 1 , for any function u = (u 1 , u 2 ) T ∈ (W 3,∞ (Ω)) 2 and any parallelogram center x e , we have (2) Assume that and h is enough small with any parallelogram center x e , we have (3) Under the assumption of (1) and (2) above, E and E h satisfy Eqs. (4) and (6), respectively. If E ∈ (W 3,∞ (Ω)) 2 and h is enough small, then we have For simplicity, we use the Ndof and I to denote the number of degrees of freedom and the identity matrix, respectively. We choose the analytic solution and compute F accordingly.
The Example 4.1 is aimed to test the convergence result stated in Theorem 3.1. Clearly, the results shown in Table 1 are consistent with Theorem 3.1.  Example 4.2 simulates the backward wave propagation phenomenon illustrated in Fig.1, which is based on the three types a-posteriori error indicators mentioned above. It can see that wave propagates as usual before it enters the metamaterial region. Once the wave enters into the metamaterial, wave propagates backward inside the metamaterial slab. And that wave propagates as usual again when leaves the metamaterial slab.
From the point of view of suitability between the grid and the numerical solution, it can be seen that the residual type is better than the others in Fig.1. From the numerical solution of the images, some parts of the horizontal interfaces between the two media do not need to be refined, but in fact, these places are still refined when using the recovery a-posteriori error indicators.
Thus according to two horizontal interfaces formed between the media, we divide the whole area into three sub-areas, and the recovery type a-posteriori error of the whole area is consisted of these sub-areas, which are shown in Fig.2. Compared with the direct recovery type a-posteriori error of the whole area in Fig.1, adaptive meshes by using the recovery type a-posteriori error of the whole area consisted of the error for these sub-areas perform better on the interfaces during the refinement process, which achieves the expected effect.
Remark 3. In Example 4.3-4.6, we also first consider the three types a posteriori error indicators mentioned above on the whole area and then consider that according to the number of interfaces n formed between the media, the entire area can be divided into n + 1 sub-areas, which each sub-area is independent and parallel computing can be considered, the recovery type a posteriori error of the whole area is consisted of the error for these sub-areas. which can be easily derived according to the theory of article [36] while the wave propagates along the y-axis. The constant m p is called "moving parameter" and the symbol " * " denotes as the symmetric part of the matrix. In our simulation, we choose a uniform mesh with h = 0.1 and use a PML with 5 cells in thickness around the rectangular domain. 8510 Ndof (for the initial mesh), 133620 Ndof (by using η r0 K l ) after 14 refinements, 139743 Ndof (by using η r1 K l ) and 132334 Ndof (by using η r2 K l ) with the same times of 12 refinements. Then, the center point x 0 = (x 0 , y 0 ) are chosen as (−1, 1.45) with m p = −2 for the domain Ω 1 , i.e., the wave source is fixed in the upper left. Obviously, we can see that the wave moves from left to right after wave passes the moving medium area shown in Figs.6 and 7.
Finally, the constant m p is taken as 2 for domain Ω 2 and −2 for domain Ω 3 with the center point x 0 = (0, 1.45), i.e., the wave source is fixed in the upper middle. Evidently, we can see that the wave moves from the middle to left and right after wave passes the moving medium area shown in Figs.8 and 9.
Overall, it can be seen that the recovery type a-posteriori error indicators perform better.  K l ) after 18 refinements, 315182 Ndof (by using η r1 K l ) and 273473 Ndof (by using η r2 K l ) with the same times of 13 refinements.    K l ) after 21 refinements, 306934 Ndof (by using η r1 K l ) and 280265 Ndof (by using η r2 K l ) with the same times of 13 refinements. where m = θ0r R2−R1 . Here, we choose the initial mesh with 167744 Ndof and take R 1 = 0.2, R 2 = 2R 1 , θ 0 = π 2 in the simulation.
In Example 4.4, we simulate the cylindrical rotation cloak. To see how wave propagates in the rotator, we plot the electric fields E in Figs.10 and 11. Clearly, it can be observed that the wave gets distorted in the metamaterial region and note that the parameter θ 0 determines the degree of distortion of the wave. Although the residual type a-posteriori error estimator achieves the rotation effect, the corresponding grid does not match the graphs of the numerical solution at all. In contrast, the recovery type a-posteriori error indicators are always good and meshes are more refined in the metamaterial region. Figure 10. Example 4.4: The first line, the second line and the third line are the real values of E 1 , the real values of E 2 and the meshes, respectively. From left to right: 344405 Ndof (by using η r0 K l ) after 3 refinements, 344620 Ndof (by using η r1 K l ) and 332141 Ndof (by using η r2 K l ) with the same times of 17 refinements.
Example 4.5. We consider the parameters 1, otherwise, Figure 11. Example 4.4: The first column, the second column and the third column are snapshots of numerical solution for the real values of E 1 , E 2 and adaptive meshes, respectively. The first line: 393423 Ndof after 23 refinements (by using η r1 K l ); The second line: 416438 Ndof after 21 refinements (by using η r2 K l ). and ε =    ε r cos 2 θ + ε θ sin 2 θ (ε r − ε θ ) sin θ cos θ * ε r sin 2 θ + ε θ cos 2 θ , a ≤ r ≤ c, where the constants are denoted, respectively, as, In our simulation, the parameters a, b and c are taken as 0.2, 0.4 and 0.6, respectively. There is 167744 Ndof in the initial mesh.
Example 4.5 simulates the electromagnetic concentration effect, which are shown in Figs.12 and 13. As in the previous example, although the residual type aposteriori error estimator achieves the concentration effect, the pattern of the corresponding grid is not matched with the graph of the numerical solution. The advantage of recovery type a-posteriori errors is prominent again. and Ω 2 (Γ = Ω 1 ∩ Ω 2 ) are shown in Fig.14. Here, we take R 1 = 0.2, R 2 = 2R 1 and the initial mesh with 20287 Ndof. We choose k = 20, P = [0, 1] t , d = [−1, 0] t and q(x) = 1. The parameters µ and ε are given as follows x ∈ Ω 1 , K l ) after 4 refinements, 506294 Ndof (by using η r1 K l ) after 23 refinements and 505234 Ndof (by using η r2 K l ) after 17 refinements. Figure 13. Example 4.5: Snapshots of numerical solution and adaptive meshes for the real values of E 2 . First two columns: 468957 Ndof (by using η r1 K l ) after 28 refinements; The last two columns: 656397 Ndof (by using η r2 K l ) after 25 refinements.
In Example 4.6, the colaking simulation is displayed in Figs.15 and 16. It is obvious that the recovery-based a-posteriori error estimation η r1 K l works well and the implementation is faster. Note that neither the permeability µ −1 nor the permittivity ε belongs to L ∞ (Ω), which is not under assumptions in Theorem 3.1. In spite of this, recovery-based a-posteriori error estimation η r2 K l still works for this example.  K l ) after 126 refinements, 323420 Ndof (by using η r1 K l ) after 10 refinements, 120272 Ndof (by using η r2 K l ) after 30 refinements. Figure 16. Example 4.6: Snapshots of numerical solution and adaptive meshes for the real values of E 2 . First two columns: 291690 Ndof (by using η r1 K l ) after 10 refinements; The last two columns: 78497 Ndof (by using η r2 K l ) after 34 refinements.
Remark 4. The relationship between the direction of wave source and the position of metamaterial plays a very important role in the actual simulation.

5.
Conclusions. In this paper, we used adaptive edge element method for the time-harmonic Maxwell's equation, which was not only based on the residual type a-posteriori estimation, but also based on the recovery type. And we presented the numerical examples to illustrate the reliability and efficiency of the proposed aposteriori error estimations. By adding the different wave source terms, the results were similar to those obtained by the time domain Maxwell's equations simulations [18,36]. Through the results of these examples, it can see that the recovery type a-posteriori error indicator deserves to be considered and recommended on the each sub-area.