APredictor-Corrector Finite ElementMethod for Time-Harmonic Maxwell’s Equations in Polygonal Domains

,e overall efficiency and accuracy of standard finite element methods may be severely reduced if the solution of the boundary value problem entails singularities. In the particular case of time-harmonic Maxwell’s equations in nonconvex polygonal domains Ω, H1-conforming nodal finite element methods may even fail to converge to the physical solution. In this paper, we present a new nodal finite element adaptation for solving time-harmonic Maxwell’s equations with perfectly conducting electric boundary condition in general polygonal domains. ,e originality of the present algorithm lies in the use of explicit extraction formulas for the coefficients of the singularities to define an iterative procedure for the improvement of the finite element solutions. A priori error estimates in the energy norm and in the L2 norm show that the new algorithm exhibits the same convergence properties as it is known for problems with regular solutions in the Sobolev space H2(Ω)2 in convex and nonconvex domains without the use of graded mesh refinements or any other modification of the bilinear form or the finite element spaces. Numerical experiments that validate the theoretical results are presented.


Introduction
e solution of time-harmonic Maxwell's equations on domains with geometric singularities exhibits some peculiar regularity properties that make its approximation significantly more difficult than other elliptic boundary value problems. In the case of two-dimensional domains Ω with corners and with a given right-hand side function in the space L 2 (Ω) 2 , it is known that the solution belongs in general to the Sobolev space H 2 (Ω) 2 only if the measure of the largest interior angle at the corners is strictly less than π/2 and that the solution belongs to the Sobolev space H 1 (Ω) 2 only if the domain Ω is convex (cf. [1][2][3][4]). In the case where the solution does not belong to the space H 2 (Ω) 2 standard, H 1 -conforming nodal finite element methods (cf. [5][6][7]) not only lose accuracy but also may even fail to converge to the physical solution, even with graded mesh refinements, when the domain Ω has reentrant corners (cf. [8,9]). As a consequence, the more expensive edge finite element discretization methods which are curl-conforming have to be used for the approximation [10][11][12][13].
It is always preferable to employ the more popular and widely used H 1 -conforming nodal finite element methods for the discretization of the Maxwell problem even in domains with corners. For this reason, some adaptive techniques have been proposed. e most popular of the available strategies are (1) the singular field method (cf. [9,[14][15][16]), (2) the orthogonal singular field method (cf. [14][15][16]), (3) the singular complement method (cf. [17]), and (4) the weighted regularization method (cf. [8]) and the local L 2 -projection method (cf. [18,19]). Refer also [20,21] for some other recent adaptive finite element approaches. e common future of all the abovementioned methods is that they are applied to time-harmonic Maxwell's equations on two-dimensional domains with reentrant corners in order to enforce the convergence of the nodal finite element approximations to the physical solution, but the accuracy of the computed solutions is still not optimal as one will expect for problems with regular solutions. It has been shown, see [4,22], that the singular field method combined with graded mesh refinements can yield optimal convergence. It is therefore still worthwhile to consider and study such adaptations of the standard nodal finite element discretization of the Maxwell problem on general polygonal domains that exhibit optimal accuracy. e main purpose of the present paper is fourfold: firstly, to propose extraction formulas for the computation of the coefficients of the singularities for time-harmonic Maxwell's equations in general polygonal domains; secondly, to present an adaptation of the standard nodal linear finite element approximation of the solution on quasiuniform meshes; thirdly, to show by means of a priori error estimates that the presented algorithm is efficient and the rate of convergence is optimal as it is known for problems with solutions u ∈ H 2 (Ω) 2 ; lastly, to present numerical experiments that validate the theoretical results.
Our algorithm makes use of the extraction formulas for the coefficients of the singularities, and it consists of starting with some initial values of the coefficients and then computing the coefficients and the regular part of the solution iteratively, hence the name "Predictor-corrector finite element method," refer [23] for the corresponding algorithm for Poisson's equation on polygonal domains. e present algorithm provides several advantages. e finite element spaces and the bilinear form are not modified, and mesh grading is not required. e stiffness matrices are sparse, well-conditioned, and symmetric and are generated only once at each level of refinement. e load vectors are modified additively during the corrector process a fixed number of times depending on the size of the largest angle at the corners. e coefficients of the singularities are computed independently with optimal accuracy. Existing nodal finite element software packages can easily be adapted to incorporate the present algorithm. e computed approximations not only converge to the physical solution, but they do so with optimal accuracy as it is known for problems with smooth solutions. It should be noted here that the idea of using extraction formulas for the coefficients of the singularity to define an iterative procedure for improving the finite element solution of the Poisson equation in domains with corners is well known, see, for example, [24,25]. However, the algorithm presented in Section 3.2 is different, and our focus is on rescuing the standard finite element approximation of Maxwell's equations in domains with corners. e algorithms presented in [24,25] relies on the fact that the singular solutions for the Poisson equation belong to the space H 1 for the error estimates, see [24]. In the case of Maxwell's equations in domains Ω ⊂ R 2 with reentrant corners, the singular solutions do not belong to the space H 1 . Hence, the algorithm will not converge in this case. Moreover, the embedded cutoff function in the splitting of the solution into a regular and a singular part in [24,25] is known to cause instability in the approximations, see, for example, [15,16].
is paper is organized as follows: In Section 2, we present the model boundary value problem, study the regularity and singularity properties of the solution, and give extraction formulas for the coefficients of the singularities. Section 3 contains a brief description of the usual finite element method and presents the new algorithm with the associated error estimates. In Section 4, we introduce several examples and study the convergence history for the approximated solutions and the coefficients of the singularities.
We assume that which implies also that div u � 0 in Ω. We have used in (1) the notations "curl" and "curl" to distinguish between the scalar and vectorial meanings of the curl operator defined by We observe that the operator u ⟼ curl curl u − ω 2 u is not elliptic. However, in order for us to apply the standard regularity theory for the solution of elliptic boundary value problems in two-dimensional domains with corners (cf. [29][30][31][32]) to the system (1), we consider subsequently the socalled regularized Maxwell's system which is elliptic (cf. [1,2,[14][15][16]28]). In view of (2) and the fact that Δu � curl curl u − gra d div u, the solution of (1) can be found by solving the following elliptic system: We will consider subsequently only problem (4) rather than (1). For the associated weak formulation, we introduce the Hilbert spaces equipped with the norms e weak solution of (4) is obtained by solving the variational problem.

Theorem 1.
Let Ω ⊂ R 2 be a bounded domain with boundary Γ. If Γ is of class C 1,1 or if Ω is a convex polygon, then the relation holds, and in these spaces, the norms ‖·‖ H N (Ω) and ‖·‖ H 0 (curl,div,Ω) are equivalent.
If Ω is a nonconvex polygon, then the space H N (Ω) is a proper closed subspace of the space H 0 (curl, div, Ω) with an infinite codimension. In the next subsection, we will present a full description of the H 2 -regularity properties of the solution u ∈ H 0 (curl, div, Ω) in polygonal domains.

Corner
Singularities. Suppose that Ω ⊂ R 2 is a simply connected polygonal domain; that is, the boundary Γ is the union of disjoint straight line segments Γ j , j � 1, . . . , J, numbered according to the positive orientation, that is, in anticlockwise direction. Let the endpoints of each Γ j be denoted by c j− 1 and c j , where we set c 0 � c J , and let the internal angle at c j be denoted by ω j , where 0 < ω j ≤ 2π and ω j ≠ π; that is, ω j is the angle between Γ j+1 and Γ j . Here ω j � 2π is used to model screens (cracks) and in which case Γ j+1 � Γ j . Let (r j , θ j ) denote local polar coordinates with the origin at c j , such that Γ j is on the line θ j � ω j and Γ j+1 is on the line θ j � 0 (see Figure 1). Accordingly, local Cartesian coordinates (x j , y j ) with origin at c j are defined by For subsequent use, let us introduce the following notations.

Theorem 2.
Let Ω ⊂ R 2 be a polygonal domain and let u ∈ H 0 (curl, div, Ω) be the solution of the variational problem (7) with the right-hand side function f ∈ L 2 (Ω) 2 . en, there exist unique real numbers c jl such that u can be split in a regular and a singular part in the form Moreover, there exists a constant C > 0 such that Remark 1. We note that if ω j > π, then the functions s j1 from (11) belong to the space H(curl, div, Ω) but not to the space H 1 (Ω) 2 . It follows immediately from (4) and (12) and the fact that the functions s jl are harmonic that w, the regular part of u, is the unique solution of the boundary value problem:

Mathematical Problems in Engineering
We observe from (14) that if the coefficients c jl were known, then the solution w ∈ H N (Ω) ∩ H 2 (Ω) 2 can be optimally approximated by means of the standard H 1 -conforming nodal finite element methods, irrespective of the singular nature of the exact solution u ∈ H 0 (curl, div, Ω) of (7). In this case, an approximation of u is then obtained using relation (12).

Coefficients of the Singularities.
We present subsequently extraction formulas for the coefficients c jl in (12) in the case where α jl ∉ N\ 1 { }, see (11). First, let us introduce the following notations.
For each vertex c j of Ω, we define a circular sector K j centered at c j and with radius R j and angle ω j , that is, where the radius Figure 2). Furthermore, we introduce a smooth cutoff function η j ∈ D(Ω) which depends only on the distance r j from c j by that is, supp(η j ) ⊂ K j . Finally, we define on K j the functions where , and the parameter ω 2 are taken from (4). (12) is given by the following formula: where f j � (f 1j , f 2j ) and K j , R j , r j , and θ j are defined as in (17), (15), and (10), respectively. Moreover, the estimate holds.
Proof. e derivation of formula (18) is given in [3,4], and so we omit it here for the sake of brevity. We prove subsequently the estimate (19). Let s − jl � (s 1 − jl , s 2 − jl ) be defined by It follows from (18) and the integration by parts formula the relations  Mathematical Problems in Engineering We have also used in (21) the fact that the functions s − jl are harmonic and η j varies only radially. We get from (21) with the help of Cauchy-Schwarz's inequality and relation (8) the estimates (N/B: here and elsewhere C > 0 always denotes a generic constant which may have different values at different points): We take note of the fact that □ Remark 2. We note that the situation where some singular exponents α jl are integers occurs only when α jl � 2. For example, if ω j � 3π/2, then α j3 � 2. us, suppose that α jl � 2. en, from (11), we get by direct computation that for any 0 < ϵ ≤ 2. is type of very weak singularity in H 2 (Ω) 2 does not significantly reduce the accuracy of the nodal linear finite element method which is our concern here. Consequently, no further adaptation is required for this type of singularity in order to achieve the optimal rate of convergence. However, if higher order polynomials are to be employed, then this type of singularity also reduces the accuracy.

A Predictor-Corrector Finite Element Method
In this section, we introduce and analyse an adaptation of the standard H 1 -conforming nodal linear finite element method for the numerical treatment of the variational problem (7) in polygonal domains.

e Finite Element Method and Error Estimates.
Suppose Ω ⊂ R 2 is a polygonal domain and let T h ≔ T { } h be a partition (triangulation) of the set Ω into disjoint triangles T such that the usual assumptions for H 1 -conforming finite element method are satisfied (cf. [5][6][7]). Here, Once a triangulation T h of the domain Ω has been defined, we introduce the finite dimensional finite element spaces as follows.
Definition 2. Let T h be a triangulation of the polygonal domain Ω. en, we define the space (7) is determined by formerly solving the Galerkin equation: where the sesquilinear form a(·, ·) and the linear form h(·) are taken from (7).

Remark 3.
If the polygonal domain Ω is nonconvex, then by eorem 1, the solution u ∈ H 0 (curl, div, Ω) of problem (7) may not belong to the Sobolev space H 1 (Ω) 2 and consequently cannot be approximated by the H 1 -conforming finite element method described above. In the next subsection, we introduce an adaptation of problem (28) that will yield a good approximation of the solution u ∈ H 0 (curl, div, Ω) even in nonconvex polygonal domains. e accuracy of the Galerkin solution u h ∈ V h (Ω) on quasiuniform meshes measured in the norm of the Sobolev spaces H l (Ω) (l � 0, 1) is given as follows (cf. [5][6][7]).  (28). If the solution u ∈ H 0 (curl, div, Ω) of the variational problem (7) satisfies additionally the regularity property u ∈ H σ+1 (Ω) 2 (0 < σ ≤ 1), then there exists a constant C > 0 independent of u and h such that It follows from eorem 4 that the accuracy of the standard finite element method using linear Lagrange finite Mathematical Problems in Engineering 5 elements is severely reduced in the case where the solution is not in the space H 2 (Ω) 2 , even if the domain is convex. e algorithm presented in the next subsection also yields optimal accuracy as known for problems with solutions u ∈ H 0 (curl, div, Ω) ∩ H 2 (Ω) 2 .

e Predictor-Corrector Finite Element Algorithm.
Suppose the solution u ∈ H 0 (curl, div, Ω) of the variational problem (7) has the form (12). en, once a finite element solution u h ∈ V h (Ω) has been computed from (28), we can then approximate straightforwardly the coefficients of the singularities c jl by where

Lemma 1. Let T h be a triangulation of the polygonal domain
Ω and let c jl and c jlh be defined by (18) and (30), respectively. en, there exists a constant C > 0 independent of h such that Proof. is follows immediately from relation (22).
2h ) of the singular part s in (12) by (ii) where the singular functions s jl are defined as in (11).
(2) Corrector and smoothing steps: For i � 1, 2, . . . , n, do the following: (i) Compute an approximation w (i) h of the regular part w in (14) by determining on the triangulation T h the Galerkin solution of the boundary value problem: (ii) For j � 1, 2, . . . , J and 0 < α jl < 2, α jl ≠ 1, l ∈ N, do the following: Compute improved approximations c (i) jlh for the coefficients c jl using the formula where ) of the singular part s in (12) by (3) Compute the final approximations u h and c jlh of the solution u ∈ H 0 (curl, div, Ω) and the coefficients c jl by setting

Remark 4
(a) e main feature of Algorithm 1 is the iterative approximation of the solution w ∈ H 2 (Ω) 2 of the boundary value problem (14) and the coefficients of the singularities c jl from (18) in step 2 of the algorithm. Since the existence of the coefficients c jl is guaranteed by eorem 3, on any triangulation T h , 0 < h < 1, of the domain Ω, the sequence of approximations c (i) jlh : i ∈ N converges irrespective of the starting value c (0) jlh , and we will show subsequently, see eorem 5, that c (i) jlh ⟶ c jl as i ⟶ ∞, h ⟶ 0. Accordingly, we have w (i) h ⟶ w as i ⟶ ∞, h ⟶ 0, see Lemma 3. Hence, on a quasiuniform family of triangulations T h : h ⟶ 0 , the finite element solutions u h defined by (38) converge to the physical solution u ∈ H 0 (curl, div, Ω) of the variational problem (7) even when Ω ⊂ R 2 is nonconvex, see eorem 5. is result is essentially due to the fact that the singular functions s jl in (12)

Error Analysis.
We analyse subsequently the accuracy and performance of Algorithm 1 by estimating the associated errors. We shall study the global discretization error in the energy norm and L 2 -norm, as well as the error of approximation of the coefficients of the singularities. (12) with α jl ≔ lπ/ω j . Let T h , 0 < h < h 0 ≤ 1, be a triangulation of the polygonal domain Ω. If 0 < α jl < 2, then there exists a constant C > 0 independent of h such that

Lemma 2. Let s jl be the singular function defined in
Proof. We have In local polar coordinates r j , θ j , see (10), we get by means of Cauchy-Schwarz's inequality the estimate where r 0 ≔ inf (r j ,θ j )∈T r j . Assertion (39) follows by combining (40) and (41). □ Lemma 3. Let u ∈ H 0 (curl, div, Ω) be the solution of the variational problem (7) with the right-hand side f ∈ L 2 (Ω) 2 . Let w ∈ H(curl, div, Ω) ∩ H 2 (Ω) 2 denote the regular part of u according to the splitting (12). Let the polygonal domain Ω be covered with a quasiuniform family of triangulations T h , the Galerkin solution of the boundary value problem (34).
Proof. Apply successively first Strang's lemma, Cea's lemma, Cauchy-Schwarz's inequality, eorem 4, and Lemma 2 to obtain the estimates where l(·) and l (i) (·) denote the linear forms associated with the weak formulations of the boundary value problems (14) and (34), respectively. In (43), Π h denotes the P 1 -Lagrange interpolation operator on the triangulation T h which is defined such that(Π h w)| T � Π T w for all T ∈ T h . Assertion (42) is then a consequence of relations (13) and (43). □ Theorem 5. Let u ∈ H 0 (curl, div, Ω) be the solution of the variational problem (7) with the right-hand side f ∈ L 2 (Ω) 2 .

Mathematical Problems in Engineering
Let the polygonal domain Ω be covered with a quasiuniform family of triangulations T h , 0 < h < h 0 ≤ 1, and let u h � u (n) h be the finite element approximation of u defined as in (38) after n iterations. If nα ≥ 1 with α defined as in (42), then there exists a constant C > 0 independent of h such that Proof. Proof. It follows from relation (38), Lemma 3, and Lemma 2 the inequalities where we have used in (45) the fact that and also the fact that, by construction, We also note that, for w ∈ H 1 (Ω) 2 , the norms ‖w‖ H 0 (curl,div,Ω) and ‖w‖ H 1 (Ω) 2 are equivalent. Assertion (44) follows by repeating the estimates for n − 1, n − 2, . . . , 1 { } and taking into consideration the fact that nα ≥ 1. □ Theorem 6. Let u ∈ H 0 (curl, div, Ω) be the solution of the variational problem (7) with the right-hand side f ∈ L 2 (Ω) 2 . Let w ∈ H(curl, div, Ω) ∩ H 2 (Ω) 2 denote the regular part of u according to the splitting (12). Let the polygonal domain Ω be covered with a quasiuniform family of triangulations J h , 0 < h < h 0 ≤ 1, and let u h � u (n) h and c jlh � c (n) jlh be the finite element approximations of the solution u and the singular coefficient c jl , respectively, defined as in (38) after n iterations. If nα ≥ 1 with α defined as in (42), then there exists a constant C > 0 independent of h such that Proof. Proof. In order to prove the estimate (47), we employ the Aubin-Nitsche lemma (cf. [5,6]). us, we have where z g ∈ H 0 (curl, div, Ω) is the unique solution of the variational problem (7) corresponding to the datum g ∈ L 2 (Ω) 2 and z h is its finite element approximation according to Algorithm 1. Assertion (47) follows by applying eorem 5 to estimate the errors of the first and second factors on the right-hand side of inequality (49). e error bound (48) follows by taking account of the inequality (see (22)) and the estimate (47).

Numerical Experiments
e main purpose of this section is to present some computational results carried out with Algorithm 1. Even though the examples considered may not represent any real world situations, the difficulties due to singularities in the solution are well reflected. We will consider examples with known and unknown solutions in convex and nonconvex domains with only one singular corner and will study the convergence history of the approximated solutions in the energy norm and the L 2 -norm as well as the convergence history of the approximated coefficients of the singularities.
All computations have been carried out using a program package in Python developed in the department by the second author of this paper and executed on a laptop computer. us, very large systems could not be handled. However, as we shall see, the results obtained are sufficient to be able to make definite conclusions about the efficiency of the proposed algorithm. e resulting linear systems in the algorithm are solved iteratively by means of preconditioned conjugate gradient method, and all the integrals are computed using a 7-point Legendre Gauss-Lobatto formula. It was observed that more quadrature points did not change the results significantly. us, we can assume that the errors due to the iterative solution of the linear systems and numerical integration are negligible in comparison with the discretization error.
In all the examples, starting from an initial triangulation J h 0 , the finer triangulations are obtained by successively dividing each triangle into four congruent ones (see Figures 3 and 4). On each triangulation J h i , an approximate solution u (n) h i is computed after n iterations (n � 2, 3) of the algorithm, see (38). If the exact solution u ∈ H 0 (curl, div, Ω) is known, then the order of convergence ] is computed from two successive approximations by On the contrary, if the exact solution is not known, then the order of convergence ] is computed from three successive approximations by Similar formulas are used for the computations of the order of convergence for the approximated coefficients of the singularities. We shall use subsequently the abbreviation H 0 ≔ H 0 (curl, div, Ω). Experiment 1. In this first example, we consider a convex polygonal domain Ω, where the interior angle at the origin is 3π/4, see Figure 3; that is, the singular exponent is α � 4/3. e right-hand side datum f is chosen such that assumption (2) is satisfied and such that the solution of the variational problem (7) with ω 2 � − 1 is given in polar coordinates (r, θ) by with p(r, θ) � φ(r)r α cos(αθ), and with constants a 0 � 1472, a 1 � − 30240, a 2 � 272160, a 3 � − 1406160, a 4 � 4592700, a 5 � − 9828378, a 6 � 1377 8100, a 7 � − 12203460, a 8 � 6200145, and a 9 � − 1377810. e truncation function φ belongs to the space C 3 (Ω) and the solution u belongs to the space H s (Ω) (0 ≤ s < 4/3).
at is, the solution u exhibits singular behaviour near the origin. By splitting (12), the singular part s of the solution is given in polar coordinates by According to the theory, see eorem 4, for this problem the standard nodal P 1 -finite element approximation converges to the physical solution and the accuracy is of the order O(h (4/3)− l ) in the norms of H l (Ω) 2 (l � 0, 1), which is not optimal. In this situation, the traditional adaptive strategies such as those presented in [8,9,[14][15][16][17] cannot be used to improve the accuracy. Tables 1 and 2 show the error estimates ‖u − u (n) h ‖ in the L 2 -norm and the energy norm after n iterations (n � 0, 1, 2, 3) of Algorithm 1, where n � 0 means no adaptation has been used. Table 3 shows the convergence history of the approximated coefficient c (n) h i . As predicted by eorems 5 and 6, the optimal rate of convergence is attained by Algorithm 1.

Experiment 2.
Here, we consider the Γ-shape domain Ω given by Figure 4. For this domain with a nonconvex corner at the origin with interior angle 3π/2, it is known that standard nodal H 1 -conforming finite element methods do not, in general, converge to the physical solution of Maxwell's equations due to the presence of two strong nonlogarithmic singularities. We show subsequently that Algorithm 1 presents an efficient adaptation. For this purpose, we choose a right-hand side datum f such that assumption (2) is satisfied and such that the solution of the variational problem (7) with ω 2 � − 1 is given in polar coordinates (r, θ) by Mathematical Problems in Engineering where p 1 (r, θ) � φ(r)r α 1 cos α 1 θ , and where the function φ is taken from (54). According to eorem 2, the singular part s of the solution is given in polar coordinates by s � (58) Table 4 shows the convergence history of the approximated coefficients c (n) 1h and c (n) 2h after n iterations (n � 2, 3). Table 5 presents estimates of the error u − u (n) h in the L 2 -norm and the energy norm. We observe that the numerical computations confirm the theoretical predictions of eorems 5 and 6. e convergence of Algorithm 1 is essentially due to the fact that the singular function s is treated semianalytically, see (37).

Experiment 3.
e goal of this experiment is to demonstrate that Algorithm 1 is efficient even for problems with combined algebraic and logarithmic singularities, see (11). For this purpose, we consider again the Γ-shape domain, see Figure 4, and choose the right-hand side datum f � (1, 1) T and ω 2 � − 1.
Although the exact solution u ∈ H 0 (curl, div, Ω) of problem (7) with this datum is not known, by eorem 2, the solution exhibits singular behaviour near all the six vertices of the domain. In fact, near the nonconvex vertex at the origin with angle 3π/2, the solution exhibits a singular behaviour of the form where the singular coefficients c i , i � 1, 2, 3, are unknown. Table 6 shows the approximated values c (n) 1h and c (n) 2h after n iterations (n � 2, 3) of the two singular coefficients c 1 and c 12 . It is obvious that c (n) 1h ⟶ 0 as h ⟶ 0. In Table 7, we present the convergence history of the approximated solutions u (n) h after n iterations (n � 2, 3) in the L 2 -norm and the energy norm. Since the exact solution is not known, we used formula (52) to compute the rates of convergence. We observe that, with this more general example, the computational results still confirm the theoretical results of eorems 5 and 3.3. In fact, the convergence rate in the energy norm is a lot higher than expected.
is, however, is an indication that the terms curl u hL 2 (Ω) 2 and div u hL 2 (Ω) 2 in the energy norm are small. We also notice that the logarithmic      Table 7: Estimates of the errors ‖u (n) h i − u (n) h i+1 ‖ L 2 (Ω) 2 and ‖u (n) singularities in the solution do not affect the rate of convergence of the predictor-corrector P 1 -finite element method, as indicated in Remark 2.

Concluding Remarks and Perspectives
We have presented extraction formulas for the coefficients of the singularities of solutions of time-harmonic Maxwell's equations in two-dimensional domains with corners. Using the explicit formulas, we proposed a predictor-corrector P 1 -conforming finite element algorithm on quasiuniform meshes for the numerical solution. We showed by means of a priori error estimates that the proposed algorithm is very efficient and it exhibits the same accuracy as it is known for problems with regular solutions u ∈ H 2 (Ω) 2 . Several numerical computations that validate the theoretical results are also presented. e present algorithm can easily be incorporated onto any existing P 1 -conforming finite element code by simply adding a suitable quadrature formula for the computation of the coefficients of the singularities and adjusting the computation of the load vector. It should be noted that calculating the coefficients of the singularities of solutions of boundary value problems is an important problem in computational mechanics (cf. [34]). us, the present algorithm solves this problem in the case of time-harmonic Maxwell's equations in two-dimensional domains.
It should be noted that the algorithm presented here can be employed efficiently with quadratic finite elements. However, the application of higher order finite elements will require that the logarithmic singularities in the solution are also taken care of.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.