# A Predictor-Corrector Finite Element Method for Time-Harmonic Maxwell’s Equations in Polygonal Domains

Authors: Boniface Nkemzi; Jake Léonard Nkeck
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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Ω, 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. 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 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.


---

## Body

## 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 L2Ω2, it is known that the solution belongs in general to the Sobolev space H2Ω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 H1Ω2 only if the domain Ω is convex (cf. [1–4]). In the case where the solution does not belong to the space H2Ω2 standard, H1-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 usedH1-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 L2-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 solutionsu∈H2Ω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. 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 H1 for the error estimates, see [24]. In the case of Maxwell’s equations in domains Ω⊂ℝ2 with reentrant corners, the singular solutions do not belong to the space H1. 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 Section2, 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Ω⊂ℝ2 with boundary Γ≔∂Ω, representing a homogeneous isotropic medium and satisfying perfectly conducting electric boundary condition, are (cf. [1, 2, 14, 16, 26–28])(1)curlcurlu−ω2u=f,inΩ,u∧n=u1n2−u2n1=0,onΓ,where ω≠0 is a complex number, n=n1,n2 denotes the outward unit normal vector on Γ, and f=f1,f2∈L2Ω2 is a given datum representing the impressed current density. We assume that(2)divf=0inΩ,which implies also that divu=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(3)curlv=∂v2∂x1−∂v1∂x2,curlφ=∂φ∂x2,−∂φ∂x1.We observe that the operatoru⟼curlcurlu−ω2u 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 Δu=curlcurlu−graddivu, the solution of (1) can be found by solving the following elliptic system:(4)−Δu−ω2u=f,inΩ,u∧n=0,onΓ,divu=0,onΓ.We will consider subsequently only problem (4) rather than (1). For the associated weak formulation, we introduce the Hilbert spaces(5)Hcurl,div,Ω≔v∈L2Ω2:curlv,divv∈L2Ω,H0curl,div,Ω≔v∈Hcurl,div,Ω:v∧n=0onΓ,HNΩ≔v∈H1Ω2:v∧n=0onΓ,equipped with the norms(6)vHcurl,div,Ω≔vL2Ω22+curlvL2Ω2+divvL2Ω21/2,vH0curl,div,Ω≔vHcurl,div,Ω,vHNΩ≔vH1Ω2.The weak solution of (4) is obtained by solving the variational problem.Findu∈H0curl,div,Ω such that(7)au,v≔curlu,curlv+divu,divv−ω2u,v=f,v≕hv,∀v∈H0curl,div,Ω,where ⋅,⋅ denotes the usual scalar product in L2Ω. If ω2≠0 is not an eigenvalue of the Dirichlet–Laplace operator in Ω, then the variational problem (7) has for any f∈L2Ω2 a unique solution that satisfies the a priori estimate (cf. [1, 2, 14–16, 28]):(8)uL2Ω22+curluL2Ω2+divuL2Ω2≤CfL2Ω22.The following preliminary regularity result for the solutionu∈H0curl,div,Ω of the variational problem (7) will be useful for the subsequent discussion on the numerical solution (cf. [33]).Theorem 1.
LetΩ⊂ℝ2 be a bounded domain with boundary Γ. If Γ is of class C1,1 or if Ω is a convex polygon, then the relation(9)HNΩ=H0curl,div,Ω,holds, and in these spaces, the norms ∥⋅∥HNΩ and ∥⋅∥H0curl,div,Ω are equivalent.IfΩ is a nonconvex polygon, then the space HNΩ is a proper closed subspace of the space H0curl,div,Ω with an infinite codimension. In the next subsection, we will present a full description of the H2-regularity properties of the solution u∈H0curl,div,Ω in polygonal domains.


### 2.2. Corner Singularities

Suppose thatΩ⊂ℝ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 cj−1 and cj, where we set c0=cJ, and let the internal angle at cj 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 rj,θj denote local polar coordinates with the origin at cj, 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 xj,yj with origin at cj are defined by(10)xj=rjcosθj,yj=rjsinθj,orrj=xj2+yj21/2,θj=sgnyjarccosxjxj2+yj21/2.Figure 1
Polygonal domain.For subsequent use, let us introduce the following notations.For each vertexcj of Ω, we set αjl≔lπ/ωj, l∈ℕ≔1,2,…, and define the functions sjlxj,yj by(11)sjl≔rjαjl−1sinαjl−1θj,rjαjl−1cosαjl−1θj,ifαjl∉ℕ.rjαjl−1lnrjsinαjl−1θj+θjcosαjl−1θj,rjαjl−1lnrjcosαjl−1θj−θjsinαjl−1θj,ifαjl∈ℕ\1,where rj and θj are defined as in (10). We have the following result (cf. [1–4]).Theorem 2.
LetΩ⊂ℝ2 be a polygonal domain and let u∈H0curl,div,Ω be the solution of the variational problem (7) with the right-hand side function f∈L2Ω2. Then, there exist unique real numbers γjl such that u can be split in a regular and a singular part in the form(12)ux,y=wx,y+sx,y≔∑j=1J∑0<αjl≤2γjlsjlxj,yj,withw∈H2Ω2.
Moreover, there exists a constantC>0 such that(13)wH2Ω2≤CfL2Ω2.Remark 1.
We note that ifωj>π, then the functions sj1 from (11) belong to the space Hcurl,div,Ω but not to the space H1Ω2.
It follows immediately from (4) and (12) and the fact that the functions sjl are harmonic that w, the regular part of u, is the unique solution of the boundary value problem:(14)−△w−ω2w=f+ω2∑j=1J∑0<αjl≤2γjlsjl,inΩ,w∧n=−∑j=1J∑0<αjl≤2γjlsjl∧n,onΓ,divw=0,onΓ.
We observe from (14) that if the coefficients γjl were known, then the solution w∈HNΩ∩H2Ω2 can be optimally approximated by means of the standard H1-conforming nodal finite element methods, irrespective of the singular nature of the exact solution u∈H0curl,div,Ω of (7). In this case, an approximation of u is then obtained using relation (12).


### 2.3. Coefficients of the Singularities

We present subsequently extraction formulas for the coefficientsγjl in (12) in the case where αjl∉ℕ\1, see (11). First, let us introduce the following notations.For each vertexcj of Ω, we define a circular sector K˜j centered at cj and with radius Rj and angle ωj, that is,(15)K˜j≔xj,yj∈Ω:0<rj<Rj,0<θj<ωj,where the radius Rj of K˜j is chosen such that K˜j∩K˜i=∅ for each j≠i (see Figure 2). Furthermore, we introduce a smooth cutoff function ηj∈DΩ¯ which depends only on the distance rj from cj by(16)ηjxj,yj=ηjrj≔1,for0≤rj≤Rj3,0≤ηjrj≤1,forRj3≤rj≤2Rj3,0,forrj≥2Rj3,that is, suppηj⊂K˜j. Finally, we define on K˜j the functions(17)f1jxj,yj≔ηjf1+ω2u1−u1Δηj−2▽ηj⋅▽u1,f2jxj,yj≔ηjf2+ω2u2−u2Δηj−2▽ηj⋅▽u2,where f=f1,f2, u=u1,u2, and the parameter ω2 are taken from (4).Figure 2
A circular sector.Theorem 3.
LetΩ⊂ℝ2 be a polygonal domain with angles ωj. Let u=u1,u2∈H0curl,div,Ω be the solution of (7) with the right-hand side function f=f1,f2∈L2Ω2. If αjl≔lπ/ωj∉ℕ, then the coefficient γjl in (12) is given by the following formula:(18)γjl≔12ωjαjl−1∫K˜jf1jsinαjl−1θj+f2jcosαjl−1θjrj1−αjl−Rj2−2αjlrjαjl−1dxdy,where fj=f1j,f2j and K˜j, Rj, rj, and θj are defined as in (17), (15), and (10), respectively. Moreover, the estimate(19)γjl≤CfL2Ω2,holds.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 s−jl=s1−jl,s2−jl be defined by(20)s−jl≔rj1−αjl−Rj2−2αjlrjαjl−1sinαjl−1θj,rj1−αjl−Rj2−2αjlrjαjl−1cosαjl−1θj.
It follows from (18) and the integration by parts formula the relations(21)γjl=12ωjαjl−1∫K˜jηjf⋅s−jl+ω2ηju⋅s−jl−u⋅s−jlΔηj−2▽ηj⋅▽u1s1−jl−2▽ηj⋅▽u2s2−jldxdy=12ωjαjl−1∫K˜jηjf⋅s−jl+ω2ηju⋅s−jl+u⋅Δηjs−jldxdy.
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):(22)γjl≤C∫K˜f⋅s−jl+u⋅s−jl+u⋅△ηjs−jldxdy≤CfL2K˜j2+uL2K˜j2≤CfL2Ω2.
We take note of the fact that(23)∫K˜js−jl2+Δηjs−jl2dxdy≤C<∞.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. Thus, suppose that αjl=2. Then, from (11), we get by direct computation that(24)sjl∈H2−ϵΩ2,for any 0<ϵ≤2. This type of very weak singularity in H2Ω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.


## 2.1. The Model Problem

The classical time-harmonic Maxwell equations in a bounded and simply connected polygonal domainΩ⊂ℝ2 with boundary Γ≔∂Ω, representing a homogeneous isotropic medium and satisfying perfectly conducting electric boundary condition, are (cf. [1, 2, 14, 16, 26–28])(1)curlcurlu−ω2u=f,inΩ,u∧n=u1n2−u2n1=0,onΓ,where ω≠0 is a complex number, n=n1,n2 denotes the outward unit normal vector on Γ, and f=f1,f2∈L2Ω2 is a given datum representing the impressed current density. We assume that(2)divf=0inΩ,which implies also that divu=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(3)curlv=∂v2∂x1−∂v1∂x2,curlφ=∂φ∂x2,−∂φ∂x1.We observe that the operatoru⟼curlcurlu−ω2u 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 Δu=curlcurlu−graddivu, the solution of (1) can be found by solving the following elliptic system:(4)−Δu−ω2u=f,inΩ,u∧n=0,onΓ,divu=0,onΓ.We will consider subsequently only problem (4) rather than (1). For the associated weak formulation, we introduce the Hilbert spaces(5)Hcurl,div,Ω≔v∈L2Ω2:curlv,divv∈L2Ω,H0curl,div,Ω≔v∈Hcurl,div,Ω:v∧n=0onΓ,HNΩ≔v∈H1Ω2:v∧n=0onΓ,equipped with the norms(6)vHcurl,div,Ω≔vL2Ω22+curlvL2Ω2+divvL2Ω21/2,vH0curl,div,Ω≔vHcurl,div,Ω,vHNΩ≔vH1Ω2.The weak solution of (4) is obtained by solving the variational problem.Findu∈H0curl,div,Ω such that(7)au,v≔curlu,curlv+divu,divv−ω2u,v=f,v≕hv,∀v∈H0curl,div,Ω,where ⋅,⋅ denotes the usual scalar product in L2Ω. If ω2≠0 is not an eigenvalue of the Dirichlet–Laplace operator in Ω, then the variational problem (7) has for any f∈L2Ω2 a unique solution that satisfies the a priori estimate (cf. [1, 2, 14–16, 28]):(8)uL2Ω22+curluL2Ω2+divuL2Ω2≤CfL2Ω22.The following preliminary regularity result for the solutionu∈H0curl,div,Ω of the variational problem (7) will be useful for the subsequent discussion on the numerical solution (cf. [33]).Theorem 1.
LetΩ⊂ℝ2 be a bounded domain with boundary Γ. If Γ is of class C1,1 or if Ω is a convex polygon, then the relation(9)HNΩ=H0curl,div,Ω,holds, and in these spaces, the norms ∥⋅∥HNΩ and ∥⋅∥H0curl,div,Ω are equivalent.IfΩ is a nonconvex polygon, then the space HNΩ is a proper closed subspace of the space H0curl,div,Ω with an infinite codimension. In the next subsection, we will present a full description of the H2-regularity properties of the solution u∈H0curl,div,Ω in polygonal domains.


## 2.2. Corner Singularities

Suppose thatΩ⊂ℝ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 cj−1 and cj, where we set c0=cJ, and let the internal angle at cj 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 rj,θj denote local polar coordinates with the origin at cj, 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 xj,yj with origin at cj are defined by(10)xj=rjcosθj,yj=rjsinθj,orrj=xj2+yj21/2,θj=sgnyjarccosxjxj2+yj21/2.Figure 1
Polygonal domain.For subsequent use, let us introduce the following notations.For each vertexcj of Ω, we set αjl≔lπ/ωj, l∈ℕ≔1,2,…, and define the functions sjlxj,yj by(11)sjl≔rjαjl−1sinαjl−1θj,rjαjl−1cosαjl−1θj,ifαjl∉ℕ.rjαjl−1lnrjsinαjl−1θj+θjcosαjl−1θj,rjαjl−1lnrjcosαjl−1θj−θjsinαjl−1θj,ifαjl∈ℕ\1,where rj and θj are defined as in (10). We have the following result (cf. [1–4]).Theorem 2.
LetΩ⊂ℝ2 be a polygonal domain and let u∈H0curl,div,Ω be the solution of the variational problem (7) with the right-hand side function f∈L2Ω2. Then, there exist unique real numbers γjl such that u can be split in a regular and a singular part in the form(12)ux,y=wx,y+sx,y≔∑j=1J∑0<αjl≤2γjlsjlxj,yj,withw∈H2Ω2.
Moreover, there exists a constantC>0 such that(13)wH2Ω2≤CfL2Ω2.Remark 1.
We note that ifωj>π, then the functions sj1 from (11) belong to the space Hcurl,div,Ω but not to the space H1Ω2.
It follows immediately from (4) and (12) and the fact that the functions sjl are harmonic that w, the regular part of u, is the unique solution of the boundary value problem:(14)−△w−ω2w=f+ω2∑j=1J∑0<αjl≤2γjlsjl,inΩ,w∧n=−∑j=1J∑0<αjl≤2γjlsjl∧n,onΓ,divw=0,onΓ.
We observe from (14) that if the coefficients γjl were known, then the solution w∈HNΩ∩H2Ω2 can be optimally approximated by means of the standard H1-conforming nodal finite element methods, irrespective of the singular nature of the exact solution u∈H0curl,div,Ω of (7). In this case, an approximation of u is then obtained using relation (12).


## 2.3. Coefficients of the Singularities

We present subsequently extraction formulas for the coefficientsγjl in (12) in the case where αjl∉ℕ\1, see (11). First, let us introduce the following notations.For each vertexcj of Ω, we define a circular sector K˜j centered at cj and with radius Rj and angle ωj, that is,(15)K˜j≔xj,yj∈Ω:0<rj<Rj,0<θj<ωj,where the radius Rj of K˜j is chosen such that K˜j∩K˜i=∅ for each j≠i (see Figure 2). Furthermore, we introduce a smooth cutoff function ηj∈DΩ¯ which depends only on the distance rj from cj by(16)ηjxj,yj=ηjrj≔1,for0≤rj≤Rj3,0≤ηjrj≤1,forRj3≤rj≤2Rj3,0,forrj≥2Rj3,that is, suppηj⊂K˜j. Finally, we define on K˜j the functions(17)f1jxj,yj≔ηjf1+ω2u1−u1Δηj−2▽ηj⋅▽u1,f2jxj,yj≔ηjf2+ω2u2−u2Δηj−2▽ηj⋅▽u2,where f=f1,f2, u=u1,u2, and the parameter ω2 are taken from (4).Figure 2
A circular sector.Theorem 3.
LetΩ⊂ℝ2 be a polygonal domain with angles ωj. Let u=u1,u2∈H0curl,div,Ω be the solution of (7) with the right-hand side function f=f1,f2∈L2Ω2. If αjl≔lπ/ωj∉ℕ, then the coefficient γjl in (12) is given by the following formula:(18)γjl≔12ωjαjl−1∫K˜jf1jsinαjl−1θj+f2jcosαjl−1θjrj1−αjl−Rj2−2αjlrjαjl−1dxdy,where fj=f1j,f2j and K˜j, Rj, rj, and θj are defined as in (17), (15), and (10), respectively. Moreover, the estimate(19)γjl≤CfL2Ω2,holds.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 s−jl=s1−jl,s2−jl be defined by(20)s−jl≔rj1−αjl−Rj2−2αjlrjαjl−1sinαjl−1θj,rj1−αjl−Rj2−2αjlrjαjl−1cosαjl−1θj.
It follows from (18) and the integration by parts formula the relations(21)γjl=12ωjαjl−1∫K˜jηjf⋅s−jl+ω2ηju⋅s−jl−u⋅s−jlΔηj−2▽ηj⋅▽u1s1−jl−2▽ηj⋅▽u2s2−jldxdy=12ωjαjl−1∫K˜jηjf⋅s−jl+ω2ηju⋅s−jl+u⋅Δηjs−jldxdy.
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):(22)γjl≤C∫K˜f⋅s−jl+u⋅s−jl+u⋅△ηjs−jldxdy≤CfL2K˜j2+uL2K˜j2≤CfL2Ω2.
We take note of the fact that(23)∫K˜js−jl2+Δηjs−jl2dxdy≤C<∞.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. Thus, suppose that αjl=2. Then, from (11), we get by direct computation that(24)sjl∈H2−ϵΩ2,for any 0<ϵ≤2. This type of very weak singularity in H2Ω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.


## 3. A Predictor-Corrector Finite Element Method

In this section, we introduce and analyse an adaptation of the standardH1-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Ω⊂ℝ2 is a polygonal domain and let Th≔Th be a partition (triangulation) of the set Ω¯ into disjoint triangles T such that the usual assumptions for H1-conforming finite element method are satisfied (cf. [5–7]). Here, h≔maxT∈ThhT, where hT=diamT. Let ρT=supdiamB:B⊂Tis a  ball.Definition 1.
A family of triangulationsTh:h⟶0 is said to be quasiuniform if there exists a constant σ, independent of h and T, such that(25)ρT≥hTσ,for  eachT∈Th.
Once a triangulationTh of the domain Ω has been defined, we introduce the finite dimensional finite element spaces as follows.Definition 2.
LetTh be a triangulation of the polygonal domain Ω. Then, we define the space(26)XhΩ≔vh=v1h,v2h∈H1Ω2:vhT∈P1T2,where P1T denotes the space of all algebraic polynomials of degree ≤1 on T. Let Mi:i=1,…,dimXhΩ be the set of nodes of the triangulation Th. Then, we define the subspace(27)VhΩ≔vh∈XhΩ:vh∧nMi=0,for  eachMi∈Γ.
Of course, wheneverMi is a vertex of Γ, we let vhMi=0 instead of vh∧nMi=0.
The finite element approximationuh∈VhΩ of the solution u∈H0curl,div,Ω of the variational problem (7) is determined by formerly solving the Galerkin equation:(28)Finduh∈VhΩ:auh,vh=hvh,for allvh∈VhΩ,where the sesquilinear form a⋅,⋅ and the linear form h⋅ are taken from (7).Remark 3.
If the polygonal domainΩ is nonconvex, then by Theorem 1, the solution u∈H0curl,div,Ω of problem (7) may not belong to the Sobolev space H1Ω2 and consequently cannot be approximated by the H1-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∈H0curl,div,Ω even in nonconvex polygonal domains.
The accuracy of the Galerkin solutionuh∈VhΩ on quasiuniform meshes measured in the norm of the Sobolev spaces HlΩ l=0,1 is given as follows (cf. [5–7]).Theorem 4.
LetTh:h⟶0 be a quasiuniform family of triangulations of the polygonal domain Ω. Suppose uh∈VhΩ is the Galerkin solution defined by (28). If the solution u∈H0curl,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(29)u−uhHlΩ2≤Chmin2−l,σ+1−luHσ+1Ω2,l=0,1.It follows from Theorem4 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 H2Ω2, even if the domain is convex. The algorithm presented in the next subsection also yields optimal accuracy as known for problems with solutions u∈H0curl,div,Ω∩H2Ω2.


### 3.2. The Predictor-Corrector Finite Element Algorithm

Suppose the solutionu∈H0curl,div,Ω of the variational problem (7) has the form (12). Then, once a finite element solution uh∈VhΩ has been computed from (28), we can then approximate straightforwardly the coefficients of the singularities γjl by(30)γjlh≔12ωjαjl−1∫K˜jf1jhsinαjl−1θj+f2jhcosαjl−1θjrj1−αjl−Rj2−2αjlrjαjl−1dxdy,where(31)f1jh≔ηjf1+ω2u1h−u1hΔηj−2▽ηj⋅▽u1h,f2jh≔ηjf2+ω2u2h−u2hΔηj−2▽ηj⋅▽u2h.Lemma 1.
LetTh be a triangulation of the polygonal domain Ω and let γjl and γjlh be defined by (18) and (30), respectively. Then, there exists a constant C>0 independent of h such that(32)γjl−γjlh≤Cu−uhL2Ω2.Proof.
This follows immediately from relation (22).Algorithm 1.
(the predictor-corrector algorithm)(1)
Predictor step:(i)
Compute from (28) an initial approximation uh0≔uh of the solution u∈H0curl,div,Ω of the variational problem (7) even when Ω is nonconvex.(ii)
Forj=1,2,…,J and 0<αjl<2, αjl≠1, l∈ℕ, do the following:Compute initial approximationsγjlh0 of the coefficients γjl from (12) using formula (30).(i)
Compute an initial approximationsh0=s1h0,s2h0 of the singular part s in (12) by(33)sh0≔∑j=1J∑0<αjl<2γjlh0sjl,(ii)
where the singular functionssjl are defined as in (11).(2)
Corrector and smoothing steps:Fori=1,2,…,n, do the following:(i)
Compute an approximationwhi of the regular part w in (14) by determining on the triangulation Th the Galerkin solution of the boundary value problem:(34)−△w−ω2w=f+ω2shi−1,inΩ,w∧n=−shi−1∧n,onΓ,divw=0,onΓ.(ii)
Forj=1,2,…,J and 0<αjl<2, αjl≠1, l∈ℕ, do the following:Compute improved approximationsγjlhi for the coefficients γjl using the formula(35)γjlhi≔γjlhi−1+12ωjαjl−1∫K˜jf1jhisinαjl−1θj+f2jhicosαjl−1θjrj1−αjl−Rj2−2αjlrjαjl−1dxdy,where(36)f1jhi=ηjf1+ω2s1hi−1+ω2w1hi−w1hiΔηj−2▽ηj⋅▽w1hi,f2jhi=ηjf2+ω2s2hi−1+ω2w2hi−w2hiΔηj−2▽ηj⋅▽w2hi.(i)
Compute an improved approximationshi=s1hi,s2hi of the singular part s in (12) by(37)shi=∑j=1J∑0<αjl<2γjlhisjl.(3)
Compute the final approximationsuh and γjlh of the solution u∈H0curl,div,Ω and the coefficients γjl by setting(38)uh=uhn≔whn+shn,γjlh≔γjlhn.Remark 4.
(a)
The main feature of Algorithm1 is the iterative approximation of the solution w∈H2Ω2 of the boundary value problem (14) and the coefficients of the singularities γjl from (18) in step 2 of the algorithm. Since the existence of the coefficients γjl is guaranteed by Theorem 3, on any triangulation Th, 0<h<1, of the domain Ω, the sequence of approximations γjlhi:i∈ℕ converges irrespective of the starting value γjlh0, and we will show subsequently, see Theorem 5, that γjlhi⟶γjl as i⟶∞, h⟶0. Accordingly, we have whi⟶w as i⟶∞, h⟶0, see Lemma 3. Hence, on a quasiuniform family of triangulations Th:h⟶0, the finite element solutions uh defined by (38) converge to the physical solution u∈H0curl,div,Ω of the variational problem (7) even when Ω⊂ℝ2 is nonconvex, see Theorem 5. This result is essentially due to the fact that the singular functions sjl in (12) are taken care of in the algorithm explicitly, see (37).(b)
The efficiency of Algorithm1 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 Th and the load vectors are corrected only by additive terms. The main additional effort lies in the computation of the coefficients γjlhi by means of an appropriate Gauss quadrature formula.(c)
We observe that formula (35) for computing γjlhi contains only the approximated regular part whi and the singular part shi−1 of the solution u∈H0curl,div,Ω.


### 3.3. Error Analysis

We analyse subsequently the accuracy and performance of Algorithm1 by estimating the associated errors. We shall study the global discretization error in the energy norm and L2-norm, as well as the error of approximation of the coefficients of the singularities.Lemma 2.
Letsjl be the singular function defined in (12) with αjl≔lπ/ωj. Let Th, 0<h<h0≤1, be a triangulation of the polygonal domain Ω. If 0<αjl<2, then there exists a constant C>0 independent of h such that(39)sjlL2Ω2≤Chαjl.Proof.
We have(40)sjlL2Ω22=∑T∈ThsjlL2T22.
In local polar coordinatesrj, θj, see (10), we get by means of Cauchy–Schwarz’s inequality the estimate(41)sjlL2T22=∫Ts1jl2+s2jl2dxdy≤C∫r0r0+hTrj2αjl−1dr≤ChT2αjl,where r0≔infrj,θj∈Trj. Assertion (39) follows by combining (40) and (41).Lemma 3.
Letu∈H0curl,div,Ω be the solution of the variational problem (7) with the right-hand side f∈L2Ω2. Let w∈Hcurl,div,Ω∩H2Ω2 denote the regular part of u according to the splitting (12). Let the polygonal domain Ω be covered with a quasiuniform family of triangulations Th, 0<h<h0≤1, and let whi∈V˜hΩ≔vh∈H1Ω2:vh∧nMi=−shi−1∧nMifor each  nodeMi∈Γ be the Galerkin solution of the boundary value problem (34). Then, there exists a constant C>0 independent of h such that(42)w−whiH1Ω2≤ChfL2Ω2+hα∑j=1J∑0<αjl<2γjl−γjlhi−1,where α≔inf0<αjl<2:j=1,…,J;l∈ℕ, see (11).Proof.
Apply successively first Strang’s lemma, Cea’s lemma, Cauchy–Schwarz’s inequality, Theorem4, and Lemma 2 to obtain the estimates(43)w−whiH1Ω2≤Cinfvh∈V˜hΩw−vhH1Ω2+supzh∈V˜hΩlzh−lizhzhH1Ω2≤Cw−ΠhwH1Ω2+∑j=1J∑0<αjl<2γjl−γjlhi−1sjlL2Ω2≤C⊸wH2Ω2+hα∑j=1J∑0<αjl<2γjl−γjlhi−1,where l⋅ and li⋅ denote the linear forms associated with the weak formulations of the boundary value problems (14) and (34), respectively. In (43),Πh denotes the P1-Lagrange interpolation operator on the triangulation Th which is defined such that ΠhwT=ΠTw for all T∈Th. Assertion (42) is then a consequence of relations (13) and (43).Theorem 5.
Letu∈H0curl,div,Ω be the solution of the variational problem (7) with the right-hand side f∈L2Ω2. Let the polygonal domain Ω be covered with a quasiuniform family of triangulations Th, 0<h<h0≤1, and let uh=uhn 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(44)u−uhH0curl,div,Ω≤ChfL2Ω2.Proof.
Proof. It follows from relation (38), Lemma 3, and Lemma 2 the inequalities(45)u−uhnH0curl,div,Ω≤Cw−whnH1Ω2+∑j=1J∑0<αjl<2γjl−γjlhnsjlH0curl,div,Ω≤hfL2Ω2+hα∑j=1J∑0<αjl<2γjl−γjlhn−1+hα∑j=1J∑0<αjl<2γjl−γjlhn≤ChfL2Ω2+hα∑j=1J∑0<αjl<2γjl−γjlhn−1≤ChfL2Ω2+hαu−uhn−1H0curl,div,Ω,where we have used in (45) the fact that(46)sjlH0curl,div,Ω=sjlL2Ω2≤Chαjl,and also the fact that, by construction, γjl−γjlhn≤γjl−γjlhn−1. We also note that, for w∈H1Ω2, the norms wH0curl,div,Ω and wH1Ω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.
Letu∈H0curl,div,Ω be the solution of the variational problem (7) with the right-hand side f∈L2Ω2. Let w∈Hcurl,div,Ω∩H2Ω2 denote the regular part of u according to the splitting (12). Let the polygonal domain Ω be covered with a quasiuniform family of triangulations Jh, 0<h<h0≤1, and let uh=uhn and γjlh=γjlhn be the finite element approximations of the solution u and the singular coefficient γ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(47)u−uhL2Ω2≤Ch2fL2Ω2,(48)γjl−γjlh≤Ch2fL2Ω2.Proof.
Proof. In order to prove the estimate (47), we employ the Aubin–Nitsche lemma (cf. [5, 6]). Thus, we have(49)u−uhL2Ω2≤Cu−uhH0curl,div,Ωsupg∈L2Ω21gL2Ω2infzh∈VhΩzg−zhH0curl,div,Ω,where zg∈H0curl,div,Ω is the unique solution of the variational problem (7) corresponding to the datum g∈L2Ω2 and zh 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))(50)γjl−γjlh≤Cu−uhL2Ω2,and the estimate (47).


## 3.1. The Finite Element Method and Error Estimates

SupposeΩ⊂ℝ2 is a polygonal domain and let Th≔Th be a partition (triangulation) of the set Ω¯ into disjoint triangles T such that the usual assumptions for H1-conforming finite element method are satisfied (cf. [5–7]). Here, h≔maxT∈ThhT, where hT=diamT. Let ρT=supdiamB:B⊂Tis a  ball.Definition 1.
A family of triangulationsTh:h⟶0 is said to be quasiuniform if there exists a constant σ, independent of h and T, such that(25)ρT≥hTσ,for  eachT∈Th.
Once a triangulationTh of the domain Ω has been defined, we introduce the finite dimensional finite element spaces as follows.Definition 2.
LetTh be a triangulation of the polygonal domain Ω. Then, we define the space(26)XhΩ≔vh=v1h,v2h∈H1Ω2:vhT∈P1T2,where P1T denotes the space of all algebraic polynomials of degree ≤1 on T. Let Mi:i=1,…,dimXhΩ be the set of nodes of the triangulation Th. Then, we define the subspace(27)VhΩ≔vh∈XhΩ:vh∧nMi=0,for  eachMi∈Γ.
Of course, wheneverMi is a vertex of Γ, we let vhMi=0 instead of vh∧nMi=0.
The finite element approximationuh∈VhΩ of the solution u∈H0curl,div,Ω of the variational problem (7) is determined by formerly solving the Galerkin equation:(28)Finduh∈VhΩ:auh,vh=hvh,for allvh∈VhΩ,where the sesquilinear form a⋅,⋅ and the linear form h⋅ are taken from (7).Remark 3.
If the polygonal domainΩ is nonconvex, then by Theorem 1, the solution u∈H0curl,div,Ω of problem (7) may not belong to the Sobolev space H1Ω2 and consequently cannot be approximated by the H1-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∈H0curl,div,Ω even in nonconvex polygonal domains.
The accuracy of the Galerkin solutionuh∈VhΩ on quasiuniform meshes measured in the norm of the Sobolev spaces HlΩ l=0,1 is given as follows (cf. [5–7]).Theorem 4.
LetTh:h⟶0 be a quasiuniform family of triangulations of the polygonal domain Ω. Suppose uh∈VhΩ is the Galerkin solution defined by (28). If the solution u∈H0curl,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(29)u−uhHlΩ2≤Chmin2−l,σ+1−luHσ+1Ω2,l=0,1.It follows from Theorem4 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 H2Ω2, even if the domain is convex. The algorithm presented in the next subsection also yields optimal accuracy as known for problems with solutions u∈H0curl,div,Ω∩H2Ω2.


## 3.2. The Predictor-Corrector Finite Element Algorithm

Suppose the solutionu∈H0curl,div,Ω of the variational problem (7) has the form (12). Then, once a finite element solution uh∈VhΩ has been computed from (28), we can then approximate straightforwardly the coefficients of the singularities γjl by(30)γjlh≔12ωjαjl−1∫K˜jf1jhsinαjl−1θj+f2jhcosαjl−1θjrj1−αjl−Rj2−2αjlrjαjl−1dxdy,where(31)f1jh≔ηjf1+ω2u1h−u1hΔηj−2▽ηj⋅▽u1h,f2jh≔ηjf2+ω2u2h−u2hΔηj−2▽ηj⋅▽u2h.Lemma 1.
LetTh be a triangulation of the polygonal domain Ω and let γjl and γjlh be defined by (18) and (30), respectively. Then, there exists a constant C>0 independent of h such that(32)γjl−γjlh≤Cu−uhL2Ω2.Proof.
This follows immediately from relation (22).Algorithm 1.
(the predictor-corrector algorithm)(1)
Predictor step:(i)
Compute from (28) an initial approximation uh0≔uh of the solution u∈H0curl,div,Ω of the variational problem (7) even when Ω is nonconvex.(ii)
Forj=1,2,…,J and 0<αjl<2, αjl≠1, l∈ℕ, do the following:Compute initial approximationsγjlh0 of the coefficients γjl from (12) using formula (30).(i)
Compute an initial approximationsh0=s1h0,s2h0 of the singular part s in (12) by(33)sh0≔∑j=1J∑0<αjl<2γjlh0sjl,(ii)
where the singular functionssjl are defined as in (11).(2)
Corrector and smoothing steps:Fori=1,2,…,n, do the following:(i)
Compute an approximationwhi of the regular part w in (14) by determining on the triangulation Th the Galerkin solution of the boundary value problem:(34)−△w−ω2w=f+ω2shi−1,inΩ,w∧n=−shi−1∧n,onΓ,divw=0,onΓ.(ii)
Forj=1,2,…,J and 0<αjl<2, αjl≠1, l∈ℕ, do the following:Compute improved approximationsγjlhi for the coefficients γjl using the formula(35)γjlhi≔γjlhi−1+12ωjαjl−1∫K˜jf1jhisinαjl−1θj+f2jhicosαjl−1θjrj1−αjl−Rj2−2αjlrjαjl−1dxdy,where(36)f1jhi=ηjf1+ω2s1hi−1+ω2w1hi−w1hiΔηj−2▽ηj⋅▽w1hi,f2jhi=ηjf2+ω2s2hi−1+ω2w2hi−w2hiΔηj−2▽ηj⋅▽w2hi.(i)
Compute an improved approximationshi=s1hi,s2hi of the singular part s in (12) by(37)shi=∑j=1J∑0<αjl<2γjlhisjl.(3)
Compute the final approximationsuh and γjlh of the solution u∈H0curl,div,Ω and the coefficients γjl by setting(38)uh=uhn≔whn+shn,γjlh≔γjlhn.Remark 4.
(a)
The main feature of Algorithm1 is the iterative approximation of the solution w∈H2Ω2 of the boundary value problem (14) and the coefficients of the singularities γjl from (18) in step 2 of the algorithm. Since the existence of the coefficients γjl is guaranteed by Theorem 3, on any triangulation Th, 0<h<1, of the domain Ω, the sequence of approximations γjlhi:i∈ℕ converges irrespective of the starting value γjlh0, and we will show subsequently, see Theorem 5, that γjlhi⟶γjl as i⟶∞, h⟶0. Accordingly, we have whi⟶w as i⟶∞, h⟶0, see Lemma 3. Hence, on a quasiuniform family of triangulations Th:h⟶0, the finite element solutions uh defined by (38) converge to the physical solution u∈H0curl,div,Ω of the variational problem (7) even when Ω⊂ℝ2 is nonconvex, see Theorem 5. This result is essentially due to the fact that the singular functions sjl in (12) are taken care of in the algorithm explicitly, see (37).(b)
The efficiency of Algorithm1 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 Th and the load vectors are corrected only by additive terms. The main additional effort lies in the computation of the coefficients γjlhi by means of an appropriate Gauss quadrature formula.(c)
We observe that formula (35) for computing γjlhi contains only the approximated regular part whi and the singular part shi−1 of the solution u∈H0curl,div,Ω.


## 3.3. Error Analysis

We analyse subsequently the accuracy and performance of Algorithm1 by estimating the associated errors. We shall study the global discretization error in the energy norm and L2-norm, as well as the error of approximation of the coefficients of the singularities.Lemma 2.
Letsjl be the singular function defined in (12) with αjl≔lπ/ωj. Let Th, 0<h<h0≤1, be a triangulation of the polygonal domain Ω. If 0<αjl<2, then there exists a constant C>0 independent of h such that(39)sjlL2Ω2≤Chαjl.Proof.
We have(40)sjlL2Ω22=∑T∈ThsjlL2T22.
In local polar coordinatesrj, θj, see (10), we get by means of Cauchy–Schwarz’s inequality the estimate(41)sjlL2T22=∫Ts1jl2+s2jl2dxdy≤C∫r0r0+hTrj2αjl−1dr≤ChT2αjl,where r0≔infrj,θj∈Trj. Assertion (39) follows by combining (40) and (41).Lemma 3.
Letu∈H0curl,div,Ω be the solution of the variational problem (7) with the right-hand side f∈L2Ω2. Let w∈Hcurl,div,Ω∩H2Ω2 denote the regular part of u according to the splitting (12). Let the polygonal domain Ω be covered with a quasiuniform family of triangulations Th, 0<h<h0≤1, and let whi∈V˜hΩ≔vh∈H1Ω2:vh∧nMi=−shi−1∧nMifor each  nodeMi∈Γ be the Galerkin solution of the boundary value problem (34). Then, there exists a constant C>0 independent of h such that(42)w−whiH1Ω2≤ChfL2Ω2+hα∑j=1J∑0<αjl<2γjl−γjlhi−1,where α≔inf0<αjl<2:j=1,…,J;l∈ℕ, see (11).Proof.
Apply successively first Strang’s lemma, Cea’s lemma, Cauchy–Schwarz’s inequality, Theorem4, and Lemma 2 to obtain the estimates(43)w−whiH1Ω2≤Cinfvh∈V˜hΩw−vhH1Ω2+supzh∈V˜hΩlzh−lizhzhH1Ω2≤Cw−ΠhwH1Ω2+∑j=1J∑0<αjl<2γjl−γjlhi−1sjlL2Ω2≤C⊸wH2Ω2+hα∑j=1J∑0<αjl<2γjl−γjlhi−1,where l⋅ and li⋅ denote the linear forms associated with the weak formulations of the boundary value problems (14) and (34), respectively. In (43),Πh denotes the P1-Lagrange interpolation operator on the triangulation Th which is defined such that ΠhwT=ΠTw for all T∈Th. Assertion (42) is then a consequence of relations (13) and (43).Theorem 5.
Letu∈H0curl,div,Ω be the solution of the variational problem (7) with the right-hand side f∈L2Ω2. Let the polygonal domain Ω be covered with a quasiuniform family of triangulations Th, 0<h<h0≤1, and let uh=uhn 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(44)u−uhH0curl,div,Ω≤ChfL2Ω2.Proof.
Proof. It follows from relation (38), Lemma 3, and Lemma 2 the inequalities(45)u−uhnH0curl,div,Ω≤Cw−whnH1Ω2+∑j=1J∑0<αjl<2γjl−γjlhnsjlH0curl,div,Ω≤hfL2Ω2+hα∑j=1J∑0<αjl<2γjl−γjlhn−1+hα∑j=1J∑0<αjl<2γjl−γjlhn≤ChfL2Ω2+hα∑j=1J∑0<αjl<2γjl−γjlhn−1≤ChfL2Ω2+hαu−uhn−1H0curl,div,Ω,where we have used in (45) the fact that(46)sjlH0curl,div,Ω=sjlL2Ω2≤Chαjl,and also the fact that, by construction, γjl−γjlhn≤γjl−γjlhn−1. We also note that, for w∈H1Ω2, the norms wH0curl,div,Ω and wH1Ω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.
Letu∈H0curl,div,Ω be the solution of the variational problem (7) with the right-hand side f∈L2Ω2. Let w∈Hcurl,div,Ω∩H2Ω2 denote the regular part of u according to the splitting (12). Let the polygonal domain Ω be covered with a quasiuniform family of triangulations Jh, 0<h<h0≤1, and let uh=uhn and γjlh=γjlhn be the finite element approximations of the solution u and the singular coefficient γ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(47)u−uhL2Ω2≤Ch2fL2Ω2,(48)γjl−γjlh≤Ch2fL2Ω2.Proof.
Proof. In order to prove the estimate (47), we employ the Aubin–Nitsche lemma (cf. [5, 6]). Thus, we have(49)u−uhL2Ω2≤Cu−uhH0curl,div,Ωsupg∈L2Ω21gL2Ω2infzh∈VhΩzg−zhH0curl,div,Ω,where zg∈H0curl,div,Ω is the unique solution of the variational problem (7) corresponding to the datum g∈L2Ω2 and zh 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))(50)γjl−γjlh≤Cu−uhL2Ω2,and the estimate (47).


## 4. Numerical Experiments

The main purpose of this section is to present some computational results carried out with Algorithm1. 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 L2-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 triangulationJh0, the finer triangulations are obtained by successively dividing each triangle into four congruent ones (see Figures 3 and 4). On each triangulation Jhi, an approximate solution uhin is computed after n iterations n=2,3 of the algorithm, see (38). If the exact solution u∈H0curl,div,Ω is known, then the order of convergence ν is computed from two successive approximations by(51)ν=log2u−uhinu−uhi+1n.Figure 3
TriangulationsJh1 and Jh2.
(a)
(b)Figure 4
TriangulationsJh1 and Jh2.
(a)
(b)On the contrary, if the exact solution is not known, then the order of convergenceν is computed from three successive approximations by(52)ν=log2uhin−uhi−1nuhi+1n−uhin.Similar formulas are used for the computations of the order of convergence for the approximated coefficients of the singularities. We shall use subsequently the abbreviationH0≔H0curl,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. The 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(53)ur,θ=−163αcurlpr,θ,with(54)pr,θ=φrrαcosαθ,φr=1,if0≤r≤13,∑i=09airi,if13≤r≤23,0,ifr≤23,and with constants a0=1472, a1=−30240, a2=272160, a3=−1406160, a4=4592700, a5=−9828378, a6=13778100, a7=−12203460, a8=6200145, and a9=−1377810.
The truncation functionφ belongs to the space C3Ω¯ and the solution u belongs to the space HsΩ 0≤s<4/3. That 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(55)s=163r1/3sinθ3163r1/3cosθ3.
According to the theory, see Theorem4, for this problem the standard nodal P1-finite element approximation converges to the physical solution and the accuracy is of the order Oh4/3−l in the norms of HlΩ2 l=0,1, 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 u−uhn in the L2-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 γhin. As predicted by Theorems 5 and 6, the optimal rate of convergence is attained by Algorithm 1.Table 1
Estimates ofu−uhnL2Ω2 and rate of convergence ν.
h
u
−
u
h
0
L
2
ν
u
−
u
h
1
L
2
ν
u
−
u
h
2
L
2
ν
u
−
u
h
3
L
2
ν
h
1
3.356533
2.724650
2.73690
2.736970
h
2
1.921907
0.80
1.806122
0.59
1.747649
0.64
1.747652
0.65
h
3
0.995036
0.95
2.644371
−0.55
0.954088
0.87
0.953778
0.87
h
4
0.430117
1.21
1.743707
0.60
0.392305
1.28
0.391981
1.28
h
5
0.156936
1.45
0.231494
2.90
0.112715
1.80
0.112475
1.80
h
6
0.073838
1.09
0.112475
1.04
0.029243
1.95
0.029243
1.94
ν
exp
1.33
2.00
2.00
2.00Table 2
Estimates ofu−uhnH0curl,div,Ω and rate of convergence ν.
h
u
−
u
h
0
H
0
ν
u
−
u
h
1
H
0
ν
u
−
u
h
2
H
0
ν
u
−
u
h
3
H
0
ν
h
1
39.722314
39.725429
39.729055
39.729074
h
2
35.040322
0.18
34.939817
0.18
34.927293
0.19
34.927315
0.19
h
3
26.169409
0.42
26.276720
0.41
26.113399
0.42
26.113429
0.42
h
4
17.291051
0.60
17.803852
0.56
17.223240
0.60
17.223256
0.60
h
5
9.273734
0.90
9.194475
0.95
9.192416
0.91
9.192419
0.91
h
6
4.780459
0.97
4.681244
0.97
4.680132
0.97
4.681244
0.97
ν
exp
0.33
1.00
1.00
1.00Table 3
Estimates ofγ−γhn and rate of convergence ν.
h
# of nodes
γ
−
γ
h
1
ν
γ
−
γ
h
2
ν
γ
−
γ
h
3
ν
h
1
12
4.91084656
4.91111188
4.91111323
h
2
35
4.43159364
0.15
4.43261815
0.15
4.43261569
0.14
h
3
117
2.48117553
0.83
2.48423552
0.84
2.48423158
0.84
h
4
425
0.90078782
1.46
0.90570958
1.45
0.90570346
1.45
h
5
1617
0.24307577
1.89
0.24307189
1.90
0.24307189
1.90
h
6
6305
0.06175547
1.98
0.06185893
1.97
0.06185887
1.97
ν
exp
2.00
2.00
2.00Experiment 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 H1-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(56)ur,θ=−83α1curlp1r,θ+22α2curlp2r,θ,withα1=23,α2=43,where(57)p1r,θ=φrrα1cosα1θ,p2r,θ=φrrα2cosα2θ,and where the function φ is taken from (54).
According to Theorem2, the singular part s of the solution is given in polar coordinates by(58)s=83r−1/3sin−θ3−22r1/3sinθ383r−1/3cos−θ3−22r1/3cosθ3.
Table4 shows the convergence history of the approximated coefficients γ1hn and γ2hn after n iterations n=2,3. Table 5 presents estimates of the error u−uhn in the L2-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 s is treated semianalytically, see (37).Table 4
Estimates ofγ1−γ1hn and γ2−γ2hn and rate of convergence ν.
h
γ
1
−
γ
1
h
2
ν
γ
1
−
γ
1
h
3
ν
γ
2
−
γ
2
h
2
ν
γ
2
−
γ
2
h
3
ν
h
1
10.08844
10.088448
0.982767
0.982952
h
2
1.524068
2.73
1.516919
2.73
0.289088
1.77
0.289307
1.76
h
3
1.402713
0.12
1.395048
0.12
0.184003
0.65
0.183986
0.65
h
4
0.543154
1.37
0.534755
1.38
0.061075
1.59
0.061071
1.59
h
5
0.148301
1.87
0.137676
1.96
0.015836
1.95
0.015833
1.95
h
6
0.043563
1.77
0.035073
1.97
0.003997
1.99
0.003995
1.99
ν
e
x
p
2.00
2.00
2.00
2.00Table 5
Estimates ofu−uhnL2Ω2 and u−uhnH0curl,div,Ω and rate of convergence ν.
h
u
−
u
h
2
L
2
ν
u
−
u
h
3
L
2
ν
u
−
u
h
2
H
0
ν
u
−
u
h
3
H
0
ν
h
1
10.829945
10.829694
87.836261
87.836415
h
2
3.571482
1.60
3.507089
1.62
70.192288
0.32
70.193940
0.32
h
3
1.980387
0.85
1.785774
0.97
44.985429
0.64
44.983218
0.64
h
4
0.907800
1.12
0.648530
1.46
27.067982
0.73
27.049539
0.73
h
5
0.422822
1.10
0.150883
2.10
15.473016
0.80
13.925060
0.95
h
6
0.152655
1.47
0.038588
1.97
7.214130
1.10
7.110957
0.97
ν
exp
2.00
2.00
1.00
1.00Experiment 3.
The goal of this experiment is to demonstrate that Algorithm1 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,1T and ω2=−1. Although the exact solution u∈H0curl,div,Ω 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 3π/2, the solution exhibits a singular behaviour of the form(59)s=γ1r−1/3sin−θ3+γ2r1/3sinθ3+γ3rlnrsinθ+θcosθγ1r−1/3cos−θ3+γ2r1/3cosθ3+γ3rlnrcosθ−θsinθ,where the singular coefficients γi, i=1,2,3, are unknown.
Table6 shows the approximated values γ1hn and γ2hn after n iterations n=2,3 of the two singular coefficients γ1 and γ12. It is obvious that γ1hn⟶0 as h⟶0. In Table 7, we present the convergence history of the approximated solutions uhn after n iterations n=2,3 in the L2-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 curluhL2Ω2 and divuhL2Ω2 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 P1-finite element method, as indicated in Remark 2.Table 6
Approximated coefficientsγ1hn and γ2hn.
h
i
γ
1
h
2
γ
1
h
3
γ
2
h
2
γ
2
h
3
h
1
−
0.5602
⋅
10
−
3
−
0.5680
⋅
10
−
3
0.451145
0.451186
h
2
−
0.2727
⋅
10
−
3
−
0.2750
⋅
10
−
3
0.445977
0.445976
h
3
0.5294
⋅
10
−
4
0.4981
⋅
10
−
4
0.445461
0.445461
h
4
−
0.9829
⋅
10
−
6
−
0.9847
⋅
10
−
6
0.445911
0.445911
h
5
0.2076
⋅
10
−
5
−
0.2050
⋅
10
−
5
0.445992
0.445992
h
6
−
0.3164
⋅
10
−
6
−
0.3184
⋅
10
−
6
0.446018
0.446018Table 7
Estimates of the errorsuhin−uhi+1nL2Ω2 and uhin−uhi+1nH0.
h
i
u
h
i
2
−
u
h
i
+
1
2
L
2
ν
u
h
i
3
−
u
h
i
+
1
3
L
2
ν
u
h
i
2
−
u
h
i
+
1
2
H
0
ν
u
h
i
3
−
u
h
i
+
1
3
H
0
ν
h
2
0.1440
⋅
10
−
1
0.1436
⋅
10
−
1
0.3707
⋅
10
−
1
0.3705
⋅
10
−
1
h
3
0.3760
⋅
10
−
2
1.94
0.3750
⋅
10
−
2
1.94
0.1019
⋅
10
−
11.86
0.1018
⋅
10
−
1
1.86
h
4
0.9463
⋅
10
−
3
1.99
0.9462
⋅
10
−
3
1.99
0.2716
⋅
10
−
2
1.91
0.2716
⋅
10
−
2
1.91
h
5
0.2460
⋅
10
−
3
1.94
0.2459
⋅
10
−
3
1.94
0.7364
⋅
10
−
3
1.88
0.7363
⋅
10
−
3
1.88
h
6
0.6240
⋅
10
−
4
1.98
0.6238
⋅
10
−
4
1.98
0.1946
⋅
10
−
3
1.92
0.1946
⋅
10
−
3
1.92
ν
exp
2.00
2.00
1.00
1.00


## 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-correctorP1-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∈H2Ω2. Several numerical computations that validate the theoretical results are also presented.The present algorithm can easily be incorporated onto any existingP1-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.



---

*Source: 3502513-2020-02-22.xml*