Integral Transformations to Handle Corner Function Singularities

Recent approaches to model corner singularities in electromagnetic analysis require a treatment substantially different from that of edge singularities. In this article, new algorithms are proposed for handling the combination of corner singularities and Green’s function singularities on quadrilateral cells in method-of-moments procedures.


I. INTRODUCTION
T HE Method of Moments (MoM) and the finite element method (FEM) often use vector basis functions defined on subdomains of simple shape [1]. It is common practice for these functions to be defined in a parent ξ -space and then mapped by coordinate transformations onto the 3-D x-y-z child-space of the observer [1]. The mapping formulas from parent to child coordinates are usually of the polynomial type so that it is simple to turn a parent into child coordinates. When the structures to be modeled contain edges or corners, it is possible to augment the set of bounded basis functions (normally vector polynomials) with an appropriate set of unbounded or singular basis functions that incorporate the singular behavior of the fields in the neighborhood of the edges [1], [2] and corners [3]. In these cases, a second nonlinear mapping from the parent to a grandparent ζ -space becomes necessary to perform the required numerical integrals [2]. In addition to these two mappings, integral equation MoM applications, such as those discussed in [3]- [6], require further transformations of variables to cancel Green's function singularity from the self and near-self integrals [7]- [25]. In the following, we propose new algorithms for transforming integrals of the most common types of corner singularities, such as those that occur on the tip of infinitely thin, not necessarily flat, plates of the type considered in [3]. The electromagnetic analysis of structures with plates typically requires the MoM solution of the electric field integral equation (EFIE), and we limit our consideration to Green's function singularities for that equation and conducting plates.
The following discussion considers a single quadrilateral cell with edges locally numbered counterclockwise from one to four. The edge order number is used as a subscript to distinguish the four parent variables ξ j , ξ j +1 , ξ j +2 , and ξ j +3 that describe the parent square cell [1] 0 ≤ ξ j +1 ≤ 1, These are dummy subscripts counted modulo 4 defined such that, in the parent space, the coordinates of the singular corner are always ξ j +1 = ξ j +2 = 0. The parent-cell is mapped onto a (curved) quadrilateral cell of the child-space. An observer fixed at some location r in the child space corresponds to a local observer parent coordinate and to a local observer grandparent coordinate In general, the u and v coordinates of the observer are obtained by solving a nonlinear system, which is not a straightforward problem apart from the very simple case of bilinear mapping from parent to child space. Basis functions able to model field singularities at the cell corner ξ j +1 = ξ j +2 = 0 are the product of a bounded function and a singular factor of the form [3] S(ξ, ν) = 1 with singularity coefficient ν restricted to the range 0 ≤ ν < 1/2 to ensure the integrability of (4) on regions that span to the corner. Various vector plots of singular basis functions of this kind are shown in [3] together with 3-D plots of their divergence. S(ξ, ν) = 1 leads us immediately to the grandparent space mapped onto the parent space by the transformation This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Equation (5) maps the parent square cell onto a convex grandparent cell; the greater the value of ν, the more the cell is distorted. Whenever the integral does not involve other singular factors besides (4), it is convenient to recast the integral in the grandparent domain and use the Jacobian of (6) dξ j +1 dξ j +2 = 4 1 − 2ν to cancel the corner function singularity using In integral equation applications, we must also deal with Green's function singularity that depends on the observer position. Because of (5), the local observer grandparent coordinates (u g , v g ) are complex when the observer parent coordinates (u, v) are outside the first quadrant, i.e., for u or v negative. Therefore, we use (5) only to deal with the case of an observer in the first quadrant (i.e., for u ≥ 0 and v ≥ 0), while, as previously discussed in [2], we use transformations other than (5) to deal with observer cases in one of the remaining three quadrants, thus working in grandparent spaces other than (5).
To introduce the reader to the subject, Section II considers the case of an observer not in the source cell and not very close to the source cell. In these cases, the MoM integral evaluation is performed without a mapping inversion.
Situations that require mapping reversal are dealt with in Sections III-VI. Sections III and IV consider the important case of an observer in the first quadrant, inside or in the neighborhood of the cell. In that situation, the singularity of the basis function is canceled by means of (5), and Green's function singularity is canceled through a pseudopolar transformation. The following sections then deal in decreasing order of importance with the less frequent cases of the observer in the other quadrants. In particular, Section V considers the cases of the observer in the fourth or second quadrant, while Section VI considers the case of the observer in the third quadrant, that is, for u < 0 and v < 0, which is a case that is expected to occur rarely in applications. Validation and numerical results are presented and discussed in Section VII.

II. OBSERVER NOT IN THE SOURCE CELL AND NOT VERY CLOSE TO THE SOURCE CELL
After mapping to the parent domain, MoM source integrals over the child-patch have the form where f (ξ ) contains the rest of the corner basis function J ξ (ξ ) is the Jacobian of the transformation from parent to child domain, and  EXTRACTED TO INTEGRATE THE BASIS FUNCTIONS  OF [3, TABLE II] AND THEIR DIVERGENCE is Green's function that depends on the child-space distance from the integration point r and the observation point r, which, at most, contains a singularity of the 1/|r − r | order.
(Recall that the Jacobian J ξ (ξ ) is a factor in the denominator of the divergence-conforming functions and also of their divergence [1].) As for the singularity factors S(ξ, ν), Table I reports, for example, the singularity factors extracted from the basis functions defined in [3, Tab. II] and the corresponding ν from (4) to use the integration schemes of this article. The basis functions in [3] are for salient, perfect electric conducting (PEC) sectors, and Table I considers only the functions with leading order singularity in the range 0 < ν o < 1. Now, by using (5)-(7), the integral (9) can be expressed in the form In the event that the observer is sufficiently far from the source cell so that Green's function singularity is not an issue, the integration in (12) can be performed by quadrature over the grandparent cell without further transformation. If ν = 0, the integration domain of (12) is a square, and the problem is simplified to the point of not having to discuss it. Otherwise, (12) is reduced to the sum of three subintegrals by noticing that the eastern edge {ξ j +1 = 1; ξ j +2 ∈ [0, 1]} of the parent cell is mapped by the curved (grandparent) segment The north edge has a similar mapping; just replace subscript 2 by subscript 1 in (13), and vice versa. However, integrating on three differently sized domains, two of which have curved edges, may be time-consuming. As an alternative, with reference to Fig. 1, we prefer to introduce (pseudo)polar coordinates (λ, σ ), centered at ξ j +1 = ξ j +2 = 0 for λ = 0, and divide the parent square into two isometric subtriangles In the event that Green's function singularity is not an issue, we integrate in polar (λ, σ ) coordinates on two subtriangles that span the parent cell. The figure shows the λ = const. and σ = const. coordinate lines (for σ = 0.2 to 0.8, step = 0.2, and λ = 0.4 to 0.9, step = 0.1) in the case of ν = 1/3.
(T A , T B ) mapped by the same square domain {0 ≤ λ ≤ 1; 0 ≤ σ ≤ 1} of the polar space by the transformation formulas with α = 1 1 − 2ν (16) so to rewrite (9) as the sum of only two subintegrals It is clear that the σ -integrals in the polar frame are delicate and for those there is the need to adequately increase the order of the quadrature formula. As a final remark, we note that MoM testing integrals have a form similar to (12), simply obtained by replacing Green's function G with the expression of the field radiated from the same or any other basis function, based on the same or on a different cell. Hence, MoM testing integrals may be evaluated by using the same technique summarized by (17)- (19). We note that it is not necessary to use the Galerkin approach for testing; it may sometimes be more convenient to use a fully polynomial set of independent testing functions instead of the singular functions.

III. SELF AND NEAR SELF SOURCE INTEGRALS FOR OBSERVER IN THE FIRST QUADRANT
When the observation point, that is Green's function singular point, lies in the first quadrant of the parent space, the nonlinear transformation (5) used to cancel the singularity of (4) brings the singularity of Green's function closer to the edge of the grandparent cell. This is not a problem if the observer is located inside the cell, but, if, instead, it is outside, it becomes more difficult to calculate the source integral because the near-self zone of the grandparent domain is mapped by a much larger region of the parent space.
Green's singularity may be canceled as in [2] by subdividing the grandparent cell into four triangles sharing a common vertex at (u g , v g ). With reference to Fig. 2, the triangles are numbered with the number that labels the edge on which they are based although, in the following, to simplify the notation, the triangle T j +1 based on the western edge ξ j +1 = 0 is indicated by the symbol T W ; similarly, T j +2 is the south triangle T S based on the southern edge ξ j +2 = 0; and T j is the north triangle T N , and T E indicates the east triangle T j −1 .
On each grandparent triangle, we then introduce a local polar reference system (λ, σ ) centered at the observation point for λ = 1 and with σ increasing from 0 to 1 in the counterclockwise sense around the vertex (u g , v g ). The coordinate line λ = 0 corresponds to one edge of the grandparent cell, that is, the segment {λ = 0; 0 ≤ σ ≤ 1} of T j maps the j th edge of the child cell, and λ = 1 is the point that maps onto r, eventually. On each subtriangle, the grandparent coordinates are (21) and the Jacobian of the polar transformation vanishes at the observer for λ = 1 and cancels Green's function singularity. The mapping formulas introduced so far are reported in Table II together with the expression of Z 1 (σ ) and Z 2 (σ ) in each subtriangle. Classical applications of this cancellation technique evaluate the integral over the entire domain of each subtriangle as points outside the quadrilateral cell always belong to two subtriangles that have opposite-sign Jacobians on these external points. That is, the integral contribution from the region outside the integration cell is automatically canceled by subtraction when using this cancellation technique. However, when dealing with corner basis functions, it is not convenient to follow the classic recipe because, if we did this, the numerical precision of the whole integral might be seriously compromised. To consider the integral contributions associated only with the subregion of a triangle belonging to the quadrilateral cell, one must overcome the last difficulty, as illustrated in Fig. 3. This figure shows, just for the north triangle, how different cases are possible depending on how much the north triangle is distorted (or warped) by the proximity of the observation point to the north side of the cell. The simplest case occurs when the north triangle is completely outside the integration domain; in this  TABLE II   MAPPINGS USED TO NUMERICALLY EVALUATE THE SOURCE INTEGRALS case, the integral on the north triangle is set to zero by default, and it is not computed. The most complicated cases occur when one or both rectilinear sides of the north triangle in the grandparent domain intersect the north edge of the grandparent cell, as shown at the right-hand side in Fig. 3. A technique for evaluating the limits of the λ-σ integration intervals on-the-fly for each subtriangle is described in the following section.

IV. ON-THE-FLY CALCULATION OF THE EXTREMES
OF THE λ-σ INTEGRATION INTERVAL The source integral on a grandparent triangular domain is made up of two nested integrals: an external one over σ plus an inner integral over λ, which is evaluated with constant σ . If the observer is in the first quadrant of the grandparent space and outside the integration cell, the constant σ -line that constitutes the domain of the inner-integral intercepts either the north or the east curved edge of the grandparent domain at a single point and for a λ-value equal to λ max less than unity. Before proceeding to calculate λ max , it is obviously necessary to know which of the two curved edges the constant σ -line intercepts, also keeping in mind that it is possible that the constant σ -line lies completely outside the integration domain for all λ values in the [0, 1] range. Thus, to integrate only on the subregion that belongs to the grandparent domain, we reduce the upper limit of the λ-integral to λ max . In turn, λ max depends on σ and is obtained by solving a system of two coupled equations that impose the equality of the Z 1 and Z 2 coordinates of the curved edge of the grandparent domain (given in Table II)   the sigma coordinate of the intersection point on the curved edge (σ N or σ E , depending on the case in question) has been obtained.

A. Algorithm to Calculate λ max for the Inner Integral
At the intersection with the north side, it turns out that the value of the azimuthal variable on the north edge satisfies the nonlinear equation whereas, at the intersection with the east side satisfies the nonlinear equation The coefficients A, B, and C of the above equations are the functions of σ given in Table III. (Recall, however, that we solve the nonlinear equation by keeping σ constant.) The x value that satisfies g(x, σ ) = 0 in the interval 0 ≤ x ≤ 1 is obtained by application of the Newton-Raphson method because the derivatives of g(x, σ ) are analytically available. Notice how the sign of the second derivative is identical to that of the coefficient A; in fact, the second derivative vanishes only at x = 0 and for x 4 = 3/(1 − 2ν).
Equations (23) and (26) share the x = 0 solution for a constant sigma line intercepting the northeast corner. The columns on the right-hand side of Table III report the value of λ max . In the particular case of ν = 0, λ max is obtained immediately with no need to solve any nonlinear equation [2].
In the numerical implementation, it is necessary to take care not to lose precision in computing the "small x" solutions (say, the x < 10 −2 ones). For example, in the region where (A +C) is small and for small x, (23) is more subject to round-off errors than its approximation This and similar problems may be solved by doubling the precision of the routines that calculate the roots of (23) and (26). However, in our tests, we did not double the precision of the routines nor did we use the series representations, such as (29), to calculate the results discussed in the last section. The solutions of (23) and (26) are obtained by using the Newton-Raphson method [26], which also needs the derivatives (24) and (27) of the said functions. This method works only if there is at most one solution in the search interval. If there are two solutions, there is a local maximum or minimum, which is found with the same algorithm searching the zeroes of (24) and (27) with the help of the second derivative (28). The zero of the first derivative is used to divide the search interval of the x solution into two parts.

B. Algorithm to Calculate the σ -Domain of the Outer Integral
For an observer in the first quadrant of the grandparent space and for each integration subtriangle, the domain of the outer integral is usually the entire interval {0 ≤ σ ≤ 1}. There are, however, notable exceptions to this rule when evaluating the subintegrals on the north and east subtriangles, as noted earlier discussing Fig. 3. In particular, for an observer outside the grandparent cell, it happens that the constant σ lines of the north triangle intercept the north side of the cell only in one or two subintervals of {0 ≤ σ ≤ 1}, while they lie completely outside the integration cell in the complementary interval (see Fig. 4). The same holds for the intersections with the eastern border of the cell for the constant σ lines of the east triangle. It is, therefore, necessary to derive the range of admissible values of σ before evaluating the inner integral on λ.
A constant σ line of the north triangle intersects the northern cell edge at λ = 0 and for σ N = σ (that is for x * = 1 − σ ). This is the trivial solution of Table III, of no interest. Instead, to integrate over λ, we need to find the second intersections (assuming they exist) and the values of sigma where they occur. The trivial (x * ) and nontrivial solutions of (23) coincide at the extreme σ t of the sigma subinterval where there are two points of intersection, and x * becomes a double root of (23) for σ = σ t . That is, both the g function and its first derivative vanish at σ = σ t . It is clear that, at the extreme, one also gets λ max (σ t ) = 0. The extreme σ t of the subinterval is found from the slopes of the lines tangent to the north edge and the slope of the constant sigma line.
For σ in the interval 0 < σ < 1, the slope of the constant σ -line depends on σ (as Z 1 and Z 2 depend on σ ) and can be positive or negative, whereas the slope of any line tangent to a point (Z 1 , Z 2 ) on the north edge is always negative. (s tan is obtained directly from the explicit expression Z 1 = Z 2 (Z −2/ν 2 − 1) 1/4 of the north edge and the parametric representation of Z 2 (σ ) in terms of σ .) The slope of these two straight-lines is the same at σ t , that is, at the point where the trivial and nontrivial solutions of (23) coincide. Therefore, by equating (30) and (31), the extreme points of the sigma subintervals are the roots of the nonlinear equation obtained by application of the Newton-Raphson method with We have already mentioned that there may be two subintervals To settle these cases, it may be helpful to find the maximum or minimum point of f by finding the σ value for which (34) vanishes; this is quite simple because the first derivative of (σ ) is The sigma subsubintervals for which the constant σ -lines of the east triangle intersect the east edge of the cell are found in a manner similar to the one discussed earlier. We leave the interested reader to carry out all the details for effective implementation of the complete algorithm.

C. Algorithm Complexity and Execution Cost
As shown in Fig. 4, to deal with the subintegral on the north (or the east) triangle, the region of the first quadrant located outside the cell is divided into five zones. Depending on the position of the observer, the σ -domain of the outer integral can be the entire interval [0, 1], a subinterval (e.g., [0, σ 1 ] or [σ 2 , 1]), or two subintervals [0, σ 1 ] and [σ 2 , 1]. If the observer is in the region of Fig. 4 marked by the star, the integral on the north (or east) subtriangle is zero by default because the remaining three subtriangles completely cover the area of the north (or east) subtriangle belonging to the grandparent cell. The boundaries of these five zones are identified in the grandparent space by the following: 1) the two curved sides of the grandparent cell; 2) the two secants passing through the corner point 1) and (ζ 1 , ζ 2 ) = (1, 0); 3) the tangents to the curved sides of the cell at the northeast corner point and at the points (ζ 1 , ζ 2 ) = (0, 1) and (ζ 1 , ζ 2 ) = (1, 0). To make a long story short, the observer region is determined by executing a short sequence of if statements. Then, one or two nonlinear equations must be solved if the observer is in a region where it is necessary to determine the sigma-domain subintervals. Finally, in the course of the numerical integration on σ , a nonlinear equation must be solved to determine λ max associated with each σ -integration point.
The nonlinear equations may be solved using the Newton-Raphson method; in our tests, the convergence was rapid and without optimizing the root search algorithm for different but near σ values (which one can certainly do to get improved performance).
At this point, it is necessary to recall that the constant sigma lines of the west and south triangles also intercept one of the two curved edges of the grandparent domain. Which of the two edges is intercepted is again determined by performing a short sequence of if statements. The σ -domain of the integrals associated with T W and T S is the whole interval [0, 1].
As far as the cost of the algorithm is concerned, the execution time of the if statements mentioned earlier is negligible compared with the time needed to solve the nonlinear equations described in Sections IV-A and IV-B. The results of the extensive tests that we made have shown that to calculate an integral over the entire parent square-domain (for an observer located outside the cell), less than 5% of the time is spent in solving all the nonlinear equations that must be solved to define the λ-σ integration domains. The remaining 95% of the time is spent calculating the integral.

V. OBSERVER IN THE FOURTH AND SECOND QUADRANT
If the observer lies in the first quadrant, the singularities of (9) are completely canceled by the transformation of variables that correctly takes into account the value of ν. Unfortunately, this is no longer possible when the observer lies outside the first quadrant. At the end of the chain of (bounded) transformations that we will introduce, our integrands will have some unbounded derivatives, and clearly, this limits the precision achievable by standard numerical quadratures.
In the fourth quadrant, we immediately cancel the singularity along the ξ j +1 = 0 line by using the transformation which defines a (new) grandparent space we use for observer in the fourth quadrant. Subsequently, we subdivide the square (grandparent) integration domain into rectilinear triangles with a common vertex at the point of observation. Integration on each subtriangle is performed by working with the usual polar reference system (λ, σ ), and this leads us to define the auxiliary function where (1 − λ) is a Jacobian factor of the polar transformation that cancels the singularity of Green's function G. The southern triangle does not contribute to the integral because it is outside the integration domain. Similarly, the integral contribution I E of the eastern triangle T E is zero for u ≥ 1. This yields A. Integral I W on the Western Triangle T W As lambda varies, and for the grandparent coordinates of the west triangle describe the straight line joining the lower right corner of the integration square with the observation point. It is easy to see that the constant sigma lines intersect the east side of the square domain at λ E for σ ∈ [0, σ W ] and the south side for σ ∈ [σ W , 1] at λ S . The expressions of σ W , λ E , and λ S are given in Table IV. Thus, the integral on the west triangle is in general given by the sum

B. Integral I N on the Northern Triangle T N
The grandparent coordinates of the north triangle are and its constant sigma lines intercept the southern side of the square integration domain for σ N ≤ σ ≤ 1 while intercepting the eastern side for 0 ≤ σ ≤ σ N (see Table IV). The singularity of the integral on the north triangle is easily eliminated by using the variable transformations of Table IV.

C. Integral I E on the Eastern Triangle T E
The grandparent coordinates of the east triangle are As noted, I E is zero for u ≥ 1; otherwise, the integral is evaluated by performing the transformations of Table IV to get

D. Observer in the Second Quadrant
For an observer in the second quadrant, the singularity along the ξ j +2 = 0 line is canceled by using a transformation "specular" to that used to treat the case of observer in the fourth quadrant We then proceed, as for the fourth quadrant, to cancel the singularity of Green's function and the remaining others with appropriate transformations of variables. It is clear that the integration formulas in the second quadrant must be specularly symmetrical with respect to the bisector (of the first and third quadrants) to those obtained in the case of the observer in the fourth quadrant. For the sake of brevity, we omit all the intermediate results that the interested reader can easily derive using the transformations indicated in Table IV (bottom).

VI. OBSERVER IN THE THIRD QUADRANT
If the observer is in the third quadrant, we subdivide the parent domain into rectilinear triangles having a common vertex on the observation point (ξ j +1 , ξ j +1 ) = (u, v) and recast the integral I given by (9) in polar coordinates. Only the T N and T E triangles contribute to I ; the contribution of T W and T S is set to zero, so one gets A second transformation involving the auxiliary functions is used to cancel the singularity along one of the two singular borders ξ j +1 = 0 or ξ j +2 = 0. (Once again, the (1 − λ) factor in F 1 and F 2 cancels the singularity of Green's function G and is part of the Jacobian of the polar transformation.) The singularity along the other singular border is canceled by a final transformation that produces the integral expressions reported in the following and in Table V.

A. Integral I N on the Northern Triangle T N
The point (λ, σ ) of the north triangle maps the point of the parent domain. At constant sigma, for joining the observation point (u, v) to the lower left corner of the parent domain. However, in the north triangle, we use (51) only for |u| < |v|, that is, for σ between zero and one. In the interval σ ∈ [σ N , 1], with σ N given by (51) for |u| < |v| or σ N = 0 for |u| ≥ |v|, the constant σ lines of the north triangle intercept the west side (ξ j +1 = 0) of the parent cell at In the complementary sigma range [0, σ N ], the constant σ lines intercept the south side (ξ j +2 = 0) at Thus, in general, the integral on the north triangle is the sum of two contributions The singularities of I N1 , I N2 are canceled by using, from the bottom up, the variable transformations reported in Table V under each integral expression, which yield Finally, the integrals are numerically performed in polar (ρ, φ) domains using, from top to bottom, the transformations of Table V to compute F 1 and F 2 . The transformations that eventually lead us to work in the polar (ρ, φ) domains are necessary because the singular corner ξ j +1 = ξ j +2 = 0 belongs to the integration domain so that both F 1 and F 2 are singular at α = 0, σ = σ N . The innermost integrals of the final expressions (57) and (58) are over ρ.

B. Integral I E on the Eastern Triangle T E
In terms of polar variables, the T E parent coordinates are The constant σ lines of T E intercept the southern side (ξ j +2 = 0) of the parent cell at in the interval σ ∈ [0, σ E ], with σ E given in Table IV. In the complementary [σ E , 1] range, the constant σ lines intercept the western side (ξ j +1 = 0) at The singularities of are canceled by using, from the bottom up, the variable transformations reported in Table V, which yield The I E1 integral for small |v| is numerically delicate because (63) contains a singularity for σ = v. (A more accurate analysis shows that the singularity of (63) is weaker than that of the first-order pole in σ = v.) The same applies to small |u| with the integral I N1 because the singularity at σ = 1 − u has not been canceled in (57). This means that the formulas for an observer in the third quadrant given in this section are useful only for sufficiently large values of |u| and |v|. As we will see in the next section, if the observer is in the third quadrant, it is more convenient to use the simplified formulas (17)- (19), which we validate precisely with the formulas reported in this section. Finally, observe that I E2 is zero and I N2 is nonzero for |u| < |v|. Similarly, I N2 is zero and I E2 is nonzero for |u| > |v|. Thus, in the end, I is at most the sum of four double integrals. Two double integrals must be computed only in the case of u = v.

VII. VALIDATION OF THE INTEGRATION ALGORITHM
In the following, our quadrature algorithms are validated by results obtained for a square child cell of the unitary side; for this situation, J ξ (ξ ) = 1. Even if we performed several tests, we report below only the results for ν = 1/5 obtained by considering two scalar functions 2ν (66) that, in the integral (9), multiply one of the following two Green's functions G(ξ, r) chosen for convenience: The value of the singularity coefficient used to perform these tests is intermediate to the values of ν needed to model PEC plates with corner angle approximately equal to 90 • . The results obtained with our algorithm using G fake and the function A(ξ ) of (65) can be compared with those obtained by running a very simple Mathematica code [27] able to calculate, for m = n = 0, the fundamental integral 1 where 2 F 1 (a, b; c; z) is the Gaussian or ordinary hypergeometric function. The integral (69) converges for ν < (1 + m + n)/2, m ≥ 0, and n ≥ 0. Substantially, in this validation, we assume a unitary Green's function, and we change the position of the observer. In this way, depending on where the observation point is located, the size and shape of the four triangles used to break down the integral vary greatly even if the exact result (69) does not depend on the position of the observation point. Then, we calculate the relative error as a function of the position of the observer and construct relative error maps, such as the one shown in Fig. 5, obtained with ν = 0.2, m = n = 0.
Typically, the highest relative error occurs only when the observer is in the immediate vicinity of the edges of the unit-side integration cell. However, observation points very close to the boundary of the parent cell, or even on the cell boundary, are practically never involved in the calculation of MoM testing integrals. The blue whiskers in the first quadrant of the parent space that appear in the figures on the left and in the center of Fig. 5 are due to lack of precision in calculating Fig. 5. Relative error obtained using our algorithm in calculating the integral (69) for m = n = 0 and ν = 1/5 by varying the observer's position within the parent-space square {−5 ≤ ξ j +1 , ξ j +2 ≤ 5}. Along the horizontal axis, ξ j +1 varies from −5 to +5, while ξ j +2 varies along the vertical axis from −5 to +5. The three color maps are in the same logarithmic scale, with values of the color bar on the right that vary from 10 −12 to 10 −6 . Relative error smaller than 10 −12 is achieved when the observer is in the white areas. The figure on the left shows the results obtained using a 32 × 32 Gaussian quadrature, while, at the center, we use 64 × 64 Gaussian quadrature and, at right, a 128 × 128 Gaussian quadrature. CPU time increases by at least a factor of 4 from one graph to the next on the right. the extremes of the σ -domains that contribute to the outer σintegrals. These errors can be reduced by increasing the order of the quadrature. Fig. 6 shows, as the observer's position (u = v) changes along the bisecting line of the parent space, the trend of the relative error (in logarithmic scale) of The corner singularity in (70) is canceled by the particular expression of the scalar function B(ξ ) so that we calculate the (supposed exact) reference result by running a simple Mathematica code where Green's function singularity is again canceled by integrating in polar coordinates on subtriangles. It is evident in Fig. 6 that, in the less expensive integration technique of Section II, (17)- (19) can be used if the observation point is at a distance from the edge of the cell equal to at least 30% of the maximum size of the cell itself, that is, in the case of the figure, for u or v greater than 1.3 or less than −0.3. The expected relative error obtained by using the integration technique of Section II is 10 −8 when using 16 Gaussian integration points in λ and σ , and it becomes of the order of 10 −12 when using 32 integration points both in λ and σ . For observation points inside the cell domain, 128-point Gauss formulas should ensure relative errors better than 10 −9 if the observer is close to the edge of the cell. Fig. 7 shows the error computing with respect to the reference result obtained by using 256 × 256 Gaussian quadrature of the self/near-self formulas of Sections III-VI. For an observer in the third quadrant,  (17)- (19) are used if the observer is sufficiently distant from the integration cell and, therefore, is in the white area. In this case, we use a Gauss quadrature that uses 32 × 32 integration points for each subtriangle of Fig. 1.
the relative error of the results obtained using the simplified formulas of Section II with 256×256 quadrature is of the order of 10 −9 . Basically, the reference chosen for the integral (72) has no more than nine or ten exact digits if the observer is in the third quadrant. At any rate, for an observer in the third quadrant (u, v < 0), the results of Fig. 7 (left) together with those of Fig. 6 show that the simplified formulas of Section II, integrated using 32 × 32 Gauss quadrature, guarantee relative errors on the order of 10 −8 . Fig. 7 (left) shows how wise it is to use a 128 × 128 Gauss quadrature in the vicinity of an edge, that is, for u and v within the [−0.1, 0.1] and [0.9, 1.1] ranges. In fact, Fig. 7 (left) does not report the results obtained for u < 1/2 using 64×64 Gauss quadrature because, for these values of u, the observer is too close to the edge of the parent cell or in the third quadrant (with v between −1/5 and 1/10).
When the observation point is within the [−0.3, −0.1], [0.1, 0.9] and [1.1, 1.3] ranges, a 64 × 64 Gaussian quadrature suffices. To summarize, the quadrature formulas to be used depend on the position of the observer in the parent space. With reference to Fig. 8, in the dark gray area in the close vicinity of the cell border, one should use the Gauss quadrature with 128 × 128 integration points. The Gaussian quadrature with 64 × 64 integration points is used in the light gray areas inside the square cell and in the buffer zone outside the dark gray area. The simpler integration formulas (17)- (19) are used when the observer is located in the white area (i.e., outside the gray areas) with Gauss quadrature formulas using 32 × 32 integration points. In this way, the integration error is always of the order of 10 −8 , or better.
As a conclusion of the previous analysis, to apply the quadrature technique illustrated in Sections III-VI, the parentto-child cell mapping must be inverted if the observer is in the vicinity of the source cell. In this regard, we must clarify that the observer's parent coordinates (u, v) obtained by mapping Fig. 9.
The aperture angle of the bottom-left corner of the child cells drawn in the figures with solid-line borders is given in the insets. These cells are obtained by four different mappings of the same parent-cell {0 ≤ ξ j +1 , ξ j +2 ≤ 1}. The figures show, in yellow (light gray), the area obtained by mapping a larger parent square {−0.2 ≤ ξ j +1 , ξ j +2 ≤ 1.2}, which, in the child-domain, surrounds the (dashed) contourline PI=1.1. For the cells with tip angle equal to 45 • and 30 • , the dark gray zone at top right is mapped by complex parent-coordinates. Quadrature formulas that require mapping reversal cannot be used if the observer lies in the dark gray zone. reversal can assume complex (not real) values if the child cell is not a parallelogram.
Furthermore, the mapping reversal associates one child point to multiple points of the parent domain. For example, the reversal of a bilinear mapping (the simplest of the various possible mappings) always offers two possible solutions for the observer parent coordinates; that is to say, one finds two different pairs of coordinates (u, v) that map the same child point.
In practical applications, this is not a problem because our integration technique requires mapping reversal only if the observer is inside or in the immediate vicinity of the cell, that is, for "parent" coordinates u and v within the interval [−, 1 + ] with, for example, = 0.3 or lower. In these cases, unless the child cell is excessively distorted, the unique solution for the (u, v) pair with u and v both real is searched for ξ j +1 and ξ j +2 in the aforementioned [−, 1 + ] range.
To recognize if the observer is inside or in the neighborhood of the child cell, one evaluates the proximity index (PI). By joining, with straight lines, the observer to the four vertices of the child cell of area S, we get a pyramid whose side faces have area S1, S2, S3, and S4. The PI is simply given by and, evidently, one gets PI = 1 whenever the observer is internal or on the border of the child cell. The closer the PI is to the unit, the closer is the observer to the cell. Fig. 9 shows the contourline PI = 1.1 for four different child cells, with increasing aspect ratio (AR) and skewness (SK) as the corner aperture angle decreases. AR is the ratio of the longest to the shortest cell's side (it should be close to unity and less than 4). The skewness is where α and β are the maximum and minimum inner angles (in degree) of the cell, respectively (SK should be less than 0.9). For the examples in Fig. 9, mapping reversal is performed only if PI is less than or equal to 1.1, and the values of u and v are found by spanning the region {− ≤ ξ j +1 , ξ j +2 ≤ 1 + }, for = 1.2. While discussing Fig. 8, we observed that should be equal to 0.3; therefore, we now have to point out that the (modified) basis functions obtained in [3, Sec. V] with the multiplicative correction technique, and their divergences, vanish along the two sides opposite the singular corner. This dampens Green's function singularity if the observer is outside the cell but near the sides opposite the singular corner, which consequently allows a reduction of the value to be used in the MoM applications discussed in [3].

VIII. CONCLUSION
This article proposes new algorithms for handling the combination of corner and Green's function singularities on quadrilateral cells. These algorithms are validated numerically by considering scalar 2-D singular functions and then used in a companion article that uses new singular basis functions to enhance the convergence of method of moments (MoM) solutions for structures containing edges and corners. The algorithms presented here can be considered a starting point since, in the authors' belief, these singular integrals are not considered elsewhere in the open literature. They are certainly susceptible to further improvements to reduce the computation time in MoM applications.