Abstract
The 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 , -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. The 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 norm show that the new algorithm exhibits the same convergence properties as it is known for problems with regular solutions in the Sobolev space 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.
1. Introduction
The 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 , it is known that the solution belongs in general to the Sobolev space only if the measure of the largest interior angle at the corners is strictly less than and that the solution belongs to the Sobolev space only if the domain is convex (cf. [1–4]). In the case where the solution does not belong to the space standard, -conforming nodal finite element methods (cf. [5–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–13].
It is always preferable to employ the more popular and widely used -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. The most popular of the available strategies are (1) the singular field method (cf. [9, 14–16]), (2) the orthogonal singular field method (cf. [14–16]), (3) the singular complement method (cf. [17]), and (4) the weighted regularization method (cf. [8]) and the local -projection method (cf. [18, 19]). Refer also [20, 21] for some other recent adaptive finite element approaches. The 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.
The 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 ; 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. The present algorithm provides several advantages. The finite element spaces and the bilinear form are not modified, and mesh grading is not required. The stiffness matrices are sparse, well-conditioned, and symmetric and are generated only once at each level of refinement. The 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. The 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. The 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. The algorithms presented in [24, 25] relies on the fact that the singular solutions for the Poisson equation belong to the space for the error estimates, see [24]. In the case of Maxwell’s equations in domains with reentrant corners, the singular solutions do not belong to the space . 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].
This 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.
2. Time-Harmonic Maxwell’s Equations in 2D
2.1. The Model Problem
The classical time-harmonic Maxwell equations in a bounded and simply connected polygonal domain with boundary , representing a homogeneous isotropic medium and satisfying perfectly conducting electric boundary condition, are (cf. [1, 2, 14, 16, 26–28])where is a complex number, denotes the outward unit normal vector on , and is a given datum representing the impressed current density. We assume thatwhich implies also that 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 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–32]) to the system (1), we consider subsequently the so-called regularized Maxwell’s system which is elliptic (cf. [1, 2, 14–16, 28]). In view of (2) and the fact that , 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 spacesequipped with the norms
The weak solution of (4) is obtained by solving the variational problem.
Find such thatwhere denotes the usual scalar product in . If is not an eigenvalue of the Dirichlet–Laplace operator in , then the variational problem (7) has for any a unique solution that satisfies the a priori estimate (cf. [1, 2, 14–16, 28]):
The following preliminary regularity result for the solution of the variational problem (7) will be useful for the subsequent discussion on the numerical solution (cf. [33]).
Theorem 1. Let be a bounded domain with boundary . If is of class or if is a convex polygon, then the relationholds, and in these spaces, the norms and are equivalent.
If is a nonconvex polygon, then the space is a proper closed subspace of the space with an infinite codimension. In the next subsection, we will present a full description of the -regularity properties of the solution in polygonal domains.
2.2. Corner Singularities
Suppose that is a simply connected polygonal domain; that is, the boundary is the union of disjoint straight line segments , , numbered according to the positive orientation, that is, in anticlockwise direction. Let the endpoints of each be denoted by and , where we set , and let the internal angle at be denoted by , where and ; that is, is the angle between and . Here is used to model screens (cracks) and in which case . Let denote local polar coordinates with the origin at , such that is on the line and is on the line (see Figure 1). Accordingly, local Cartesian coordinates with origin at are defined by
For subsequent use, let us introduce the following notations.
For each vertex of , we set , , and define the functions bywhere and are defined as in (10). We have the following result (cf. [1–4]).
Theorem 2. Let be a polygonal domain and let be the solution of the variational problem (7) with the right-hand side function . Then, there exist unique real numbers such that can be split in a regular and a singular part in the formMoreover, there exists a constant such that
Remark 1. We note that if , then the functions from (11) belong to the space but not to the space .
It follows immediately from (4) and (12) and the fact that the functions are harmonic that , the regular part of , is the unique solution of the boundary value problem:We observe from (14) that if the coefficients were known, then the solution can be optimally approximated by means of the standard -conforming nodal finite element methods, irrespective of the singular nature of the exact solution of (7). In this case, an approximation of is then obtained using relation (12).
2.3. Coefficients of the Singularities
We present subsequently extraction formulas for the coefficients in (12) in the case where , see (11). First, let us introduce the following notations.
For each vertex of , we define a circular sector centered at and with radius and angle , that is,where the radius of is chosen such that for each (see Figure 2). Furthermore, we introduce a smooth cutoff function which depends only on the distance from bythat is, supp. Finally, we define on the functionswhere , , and the parameter are taken from (4).
Theorem 3. Let be a polygonal domain with angles . Let be the solution of (7) with the right-hand side function . If , then the coefficient in (12) is given by the following formula:where and , , , and are defined as in (17), (15), and (10), respectively. Moreover, the estimateholds.
Proof. The 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 be defined byIt follows from (18) and the integration by parts formula the relationsWe have also used in (21) the fact that the functions are harmonic and 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 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 are integers occurs only when . For example, if , then . Thus, suppose that . Then, from (11), we get by direct computation thatfor any . This type of very weak singularity in 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.
3. A Predictor-Corrector Finite Element Method
In this section, we introduce and analyse an adaptation of the standard -conforming nodal linear finite element method for the numerical treatment of the variational problem (7) in polygonal domains.
3.1. The Finite Element Method and Error Estimates
Suppose is a polygonal domain and let be a partition (triangulation) of the set into disjoint triangles T such that the usual assumptions for -conforming finite element method are satisfied (cf. [5–7]). Here, , where . Let .
Definition 1. A family of triangulations is said to be quasiuniform if there exists a constant σ, independent of h and T, such thatOnce a triangulation of the domain has been defined, we introduce the finite dimensional finite element spaces as follows.
Definition 2. Let be a triangulation of the polygonal domain . Then, we define the spacewhere denotes the space of all algebraic polynomials of degree on T. Let be the set of nodes of the triangulation . Then, we define the subspaceOf course, whenever is a vertex of , we let instead of .
The finite element approximation of the solution of the variational problem (7) is determined by formerly solving the Galerkin equation:where the sesquilinear form and the linear form are taken from (7).
Remark 3. If the polygonal domain is nonconvex, then by Theorem 1, the solution of problem (7) may not belong to the Sobolev space and consequently cannot be approximated by the -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 even in nonconvex polygonal domains.
The accuracy of the Galerkin solution on quasiuniform meshes measured in the norm of the Sobolev spaces is given as follows (cf. [5–7]).
Theorem 4. Let be a quasiuniform family of triangulations of the polygonal domain . Suppose is the Galerkin solution defined by (28). If the solution of the variational problem (7) satisfies additionally the regularity property , then there exists a constant independent of and h such that
It follows from Theorem 4 that the accuracy of the standard finite element method using linear Lagrange finite elements is severely reduced in the case where the solution is not in the space , even if the domain is convex. The algorithm presented in the next subsection also yields optimal accuracy as known for problems with solutions .
3.2. The Predictor-Corrector Finite Element Algorithm
Suppose the solution of the variational problem (7) has the form (12). Then, once a finite element solution has been computed from (28), we can then approximate straightforwardly the coefficients of the singularities bywhere
Lemma 1. Let be a triangulation of the polygonal domain and let and be defined by (18) and (30), respectively. Then, there exists a constant independent of h such that
Proof. This follows immediately from relation (22).
Algorithm 1. (the predictor-corrector algorithm) (1)Predictor step:(i)Compute from (28) an initial approximation of the solution of the variational problem (7) even when is nonconvex.(ii)For and , , , do the following: Compute initial approximations of the coefficients from (12) using formula (30).(i)Compute an initial approximation of the singular part in (12) by(ii)where the singular functions are defined as in (11).(2)Corrector and smoothing steps: For , do the following:(i)Compute an approximation of the regular part in (14) by determining on the triangulation the Galerkin solution of the boundary value problem:(ii)For and , , , do the following: Compute improved approximations for the coefficients using the formula where(i)Compute an improved approximation of the singular part in (12) by(3)Compute the final approximations and of the solution and the coefficients by setting
Remark 4. (a)The main feature of Algorithm 1 is the iterative approximation of the solution of the boundary value problem (14) and the coefficients of the singularities from (18) in step 2 of the algorithm. Since the existence of the coefficients is guaranteed by Theorem 3, on any triangulation , , of the domain , the sequence of approximations converges irrespective of the starting value , and we will show subsequently, see Theorem 5, that as , . Accordingly, we have as , , see Lemma 3. Hence, on a quasiuniform family of triangulations , the finite element solutions defined by (38) converge to the physical solution of the variational problem (7) even when is nonconvex, see Theorem 5. This result is essentially due to the fact that the singular functions in (12) are taken care of in the algorithm explicitly, see (37).(b)The efficiency of Algorithm 1 lies in the fact that the global stiffness matrix associated to the Galerkin solution of the problems in (34) is generated only once on each triangulation and the load vectors are corrected only by additive terms. The main additional effort lies in the computation of the coefficients by means of an appropriate Gauss quadrature formula.(c)We observe that formula (35) for computing contains only the approximated regular part and the singular part of the solution .
3.3. 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 -norm, as well as the error of approximation of the coefficients of the singularities.
Lemma 2. Let be the singular function defined in (12) with . Let , , be a triangulation of the polygonal domain . If , then there exists a constant independent of h such that
Proof. We haveIn local polar coordinates , , see (10), we get by means of Cauchy–Schwarz’s inequality the estimatewhere . Assertion (39) follows by combining (40) and (41).
Lemma 3. Let be the solution of the variational problem (7) with the right-hand side . Let denote the regular part of according to the splitting (12). Let the polygonal domain be covered with a quasiuniform family of triangulations , , and let be the Galerkin solution of the boundary value problem (34). Then, there exists a constant independent of h such thatwhere , see (11).
Proof. Apply successively first Strang’s lemma, Cea’s lemma, Cauchy–Schwarz’s inequality, Theorem 4, and Lemma 2 to obtain the estimateswhere and denote the linear forms associated with the weak formulations of the boundary value problems (14) and (34), respectively. In (43) denotes the -Lagrange interpolation operator on the triangulation which is defined such that for all . Assertion (42) is then a consequence of relations (13) and (43).
Theorem 5. Let be the solution of the variational problem (7) with the right-hand side . Let the polygonal domain be covered with a quasiuniform family of triangulations , , and let be the finite element approximation of defined as in (38) after n iterations. If with α defined as in (42), then there exists a constant independent of h such that
Proof. Proof. It follows from relation (38), Lemma 3, and Lemma 2 the inequalitieswhere we have used in (45) the fact thatand also the fact that, by construction, . We also note that, for , the norms and are equivalent. Assertion (44) follows by repeating the estimates for and taking into consideration the fact that .
Theorem 6. Let be the solution of the variational problem (7) with the right-hand side . Let denote the regular part of according to the splitting (12). Let the polygonal domain be covered with a quasiuniform family of triangulations , , and let and be the finite element approximations of the solution and the singular coefficient , respectively, defined as in (38) after n iterations. If with α defined as in (42), then there exists a constant independent of h such that
Proof. Proof. In order to prove the estimate (47), we employ the Aubin–Nitsche lemma (cf. [5, 6]). Thus, we havewhere is the unique solution of the variational problem (7) corresponding to the datum and is its finite element approximation according to Algorithm 1. Assertion (47) follows by applying Theorem 5 to estimate the errors of the first and second factors on the right-hand side of inequality (49). The error bound (48) follows by taking account of the inequality (see (22))and the estimate (47).
4. Numerical Experiments
The 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 -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. Thus, 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. The 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. Thus, 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 , the finer triangulations are obtained by successively dividing each triangle into four congruent ones (see Figures 3 and 4). On each triangulation , an approximate solution is computed after n iterations of the algorithm, see (38). If the exact solution is known, then the order of convergence ν is computed from two successive approximations by
(a)
(b)
(a)
(b)
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 .
Experiment 1. In this first example, we consider a convex polygonal domain , where the interior angle at the origin is , see Figure 3; that is, the singular exponent is . The right-hand side datum is chosen such that assumption (2) is satisfied and such that the solution of the variational problem (7) with is given in polar coordinates bywithand with constants , , , , , , , , , and .
The truncation function φ belongs to the space and the solution belongs to the space . That is, the solution exhibits singular behaviour near the origin. By splitting (12), the singular part of the solution is given in polar coordinates byAccording to the theory, see Theorem 4, for this problem the standard nodal -finite element approximation converges to the physical solution and the accuracy is of the order in the norms of , which is not optimal. In this situation, the traditional adaptive strategies such as those presented in [8, 9, 14–17] cannot be used to improve the accuracy. Tables 1 and 2 show the error estimates in the -norm and the energy norm after n iterations of Algorithm 1, where means no adaptation has been used. Table 3 shows the convergence history of the approximated coefficient . As predicted by Theorems 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 , it is known that standard nodal -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 such that assumption (2) is satisfied and such that the solution of the variational problem (7) with is given in polar coordinates bywhereand where the function φ is taken from (54).
According to Theorem 2, the singular part of the solution is given in polar coordinates byTable 4 shows the convergence history of the approximated coefficients and after n iterations . Table 5 presents estimates of the error in the -norm and the energy norm. We observe that the numerical computations confirm the theoretical predictions of Theorems 5 and 6. The convergence of Algorithm 1 is essentially due to the fact that the singular function is treated semianalytically, see (37).
Experiment 3. The 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 and . Although the exact solution of problem (7) with this datum is not known, by Theorem 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 , the solution exhibits a singular behaviour of the formwhere the singular coefficients , , are unknown.
Table 6 shows the approximated values and after n iterations of the two singular coefficients and . It is obvious that as . In Table 7, we present the convergence history of the approximated solutions after n iterations in the -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 Theorems 5 and 3.3. In fact, the convergence rate in the energy norm is a lot higher than expected. This, however, is an indication that the terms and in the energy norm are small. We also notice that the logarithmic singularities in the solution do not affect the rate of convergence of the predictor-corrector -finite element method, as indicated in Remark 2.
5. 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 -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 . Several numerical computations that validate the theoretical results are also presented.
The present algorithm can easily be incorporated onto any existing -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]). Thus, 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
The authors declare that they have no conflicts of interest.