Evaluation of Static Potential Integrals on Triangular Domains

Static potential integrals for constant and linear sources on triangles are derived in a straightforward way. The new representations, as presented, are robust with respect to machine evaluation in important limiting cases. The potential integrals comprise up to six functions, each dependent on the relative position and orientation (with respect to an observation point) of a vertex and edge, respectively, of the source triangle. Gradients of the potentials are derived by differentiation, thus preserving relations between the representations. Each such vertex function reveals any anomalous functional behavior near its associated vertex or edge, which is useful information for devising test integral schemes. Potential plots in the source plane of an equilateral triangle illustrate such behavior, as do similar plots for each vertex function and gradient components near their associated edge and vertex.


I. INTRODUCTION
This paper revises standard representations of static potentials and their gradients with a focus on their accurate and efficient evaluation, with special attention to critical limiting cases, for sources on planar elements with constant or linearlyvarying basis functions. These integrals have historically found use especially in singularity subtraction approaches for evaluating matrix elements in the Method of Moments (MoM) [1], [2]. We specifically illustrate the approach for triangular elements and RWG bases, but the concepts apply equally well to other (planar) elements and linearly-varying bases. They are also easily adapted to deal with curvilinear elements when approximating tangent elements are used to handle kernel singularities [3]. Our goal is to evaluate the following integral operators that appear in the electric and magnetic field integral equations (EFIE and MFIE), called L and K operators in [4], respectively, that form or serve as prototypes for most integral operators and formulations of The associate editor coordinating the review of this manuscript and approving it for publication was Mengmeng Li. interest [5]: where R is the separation between an observation point r and source point r on a triangular element T ; p (r ) is a linear (RWG) vector basis [6], [7] associated with a vertex or edge of T . The three integrals in (1) are readily identified as the usual scalar potential, vector potential, and curl of vector potential, respectively, and we refer to them as such. Note that since all source designations and material parameters are removed from (1), they may be interpreted as field potentials due to currents or charges of either electric or magnetic variety. For r on T , the integrands of (1) are singular when R = 0, i.e., when r = r; when r is near T (i.e., R is small, but non-vanishing), the integrands are said to be near-singular. The two usual approaches for evaluating such integrals are the singularity subtraction or singularity cancellation methods. In singularity subtraction, one first identifies a simplified integrand form asymptotic to the integrand near its singularities. Subtracting the term from the integrand yields a difference integrand less singular than the original (i.e., regularized or smoothed) and therefore more amenable to numerical integration. The subtracted term should be analytically integrable-or at least easily evaluated numerically. Since the analytically evaluated contribution contains the singularity, it is usually the dominant contribution to the total integral; in that case, reduced accuracy in the numerically evaluated difference integral is then often acceptable. Most often, the asymptotic form chosen is simply the static limit (k → 0) of the integrand; this choice usually suffices to accelerate numerical quadrature as long as kR = k r − r π/2. The static asymptotic forms of (1) generally contain the dominant contributions to the integrals, and usually closely approximate their real parts. Thus, assuming the bases p (r ) are (linear) RWG bases defined on a triangle T , the integrals can always be constructed in terms of the following static integrals with both constant and linear source densities [1], [2], [6]: We note there also exist higher-order singularity subtraction approaches in which the phase factor exp(−jkR) is expanded in a truncated power series in R and higher order (vector or scalar) basis functions, if present, are expressed as polynomials on T . The resulting power series approximations, when extracted from the integrand, further regularize and reduce the magnitude of the difference integral contribution; the resulting difference integral must, of course, be corrected by adding back the integrals of the extracted series terms, preferably evaluated in closed form as described in [6], [8], [9]. Regardless of how many leading power series terms in the phase factor are removed, however, the resulting difference integrand often remains non-analytic in the sense that at least some higher order derivative of the integrand is singular, thus eventually limiting the convergence rate of (unweighted) Gauss quadrature rules. Usually, the removed higher order integrals are evaluated using recursion formulae [8], for which the expressions presented here for the integrals (2) can be used for initialization. Unfortunately, the expressions often involve indeterminate limits and bounded or unbounded function or derivative singularities that appear not only at sources or their boundaries, but also as canceling or removable singularities far removed from source regions. Such anomalies often create numerical difficulties that either spoil the accuracy or interrupt the flow of MoM matrix entry calculations, with the problem becoming more critical as the number of quadrature points increases to obtain higher accuracy, and observation points more closely approach critical points or regions. In contrast to singularity subtraction methods, singularity cancellation methods [10]- [14] do not require the existence of closed forms for the integrals approximating integrands of (1); instead they introduce variable transforms whose Jacobians cancel any singularities-or at least regularize the integral by partially canceling rapid variations of the integrand. Unlike subtraction methods, cancellation methods often produce analytic integrands (i.e., having derivatives of all orders) that are smooth (i.e., polynomial-like). In those cases, they then provide asymptotically exponential convergence using standard (unweighted) Gauss quadrature rules. Despite the advantages of cancellation methods, they are generally more complex to implement than subtraction methods. And recently, for highly singular integrals, a hybrid approach combining subtraction of a leading singularity with singularity cancellation applied to the difference integral has proved useful [15].
The increasing need for numerically robust, closed-form static potential integral expressions with properly handled limiting conditions has been our inspiration for revisiting and revising classic results such as [1] and [2] with an eye towards providing a more robust evaluation and better understanding of these building blocks for MoM matrix element evaluation. We report the results of that study here. Section II, briefly reviews the well-known triangle subdivision strategy that guarantees dominant (near-) singularities always appear at a subtriangle vertex, allowing the same analysis to be applied to each. Potentials for both constant and linear subtriangle sources are given in closed form, and considerations regarding their evaluation and physical interpretation are presented. The concept of fundamental vertex functions is introduced. The vertex functions are element-independent functions describing the (possibly singular) behavior of potentials or their gradients near a given vertex or edge for constant or linear sources on a planar element. In Section III, rather than following the usual approach of independently evaluating potential integrals involving gradient kernels, the gradient potentials are obtained directly from the potentials of Section II by differentiation. This approach not only provides a more direct path to their evaluation, but also more closely connects their representations. Section IV summarizes the synthesis of the desired static scalar and vector potential integrals (2) from the subtriangle potentials and gradients of the earlier Sections. Section V presents and discusses numerical results beginning with plots of potentials and vector potential components near an equilateral triangle. Vertex function plots are then presented to illustrate how these static building blocks illuminate singularity behavior at triangle vertices and edges. Conclusions are presented in Section VI. Figure 1 shows an observation point r near a source triangle T with unit normaln (the caret denotes a unit vector) and with vertices and edges ordered with respect ton via the right hand rule. We assume the edge vector of the qth edge of the parent triangle is given by q = r q−1 − r q+1 (with index arithmetic performed modulo 3), the unit normal to the qth edge in the plane of the triangle isû , and r q is the position vector of the qth vertex of T . The normal projection of the observation VOLUME 8, 2020 FIGURE 1. An observation point r is projected onto the plane of triangle T , creating three subtriangles T q , q = 1, 2, 3. The projection point r 0 is labeled vertex 1 for each subtriangle. Edge q of T is a shared edge of T q (which remains as edge 1 with respect to T q ), and vertices q + 1 and q − 1 of T become vertices 2 and 3 with respect to T q . The vertex ordering around the boundaries of T and T q induces unit surface normal directions for both parent triangles and subtriangles via the right hand rule. Subtriangles share the parent's normaln if r 0 falls inside the parent triangle as shown in (a); otherwise, as in (b), one or even two of the subtriangles will lie entirely outside the parent triangle and have opposite normals (i.e., −n). point onto the plane of the source triangle is r 0 = r − dn, where d =n · (r − r q ) , q = 1, 2, or 3, is the distance of the observation point above (d > 0) or below (d < 0) the source plane. Point r 0 is then labeled vertex 1 for each of the three subtriangles, each of whose remaining vertices are endpoints of parent triangle edge q and are ordered as shown in Fig. 1. This subdivision guarantees that the static kernel's peak at r 0 always appears at vertex 1 for each subtriangle, and hence that the same scheme may be repeated for each. We also observe that, once the potential contributions of all three subtriangles are combined, the result must be independent of r 0 .

II. SUBTRIANGLE POTENTIALS FOR CONSTANT AND LINEAR SOURCE DENSITIES
For each subtriangle T q , we introduce a local rectangular coordinate system (u, , d) with origin at r 0 , associated unit vectorsû = u |u|,ˆ = | |, and height h, as shown in Figs. 1 and 2. We add the superscript index q to any of these quantities to refer to a specific subtriangle. In the remainder of this section, however, we deal with only a single subtriangle at a time, and hence, for simplicity, we temporarily drop all such superscripts. Vertices 1 , 2 , and 3 of T q are labeled such that the common edge (edge 1 of T q and edge q of T ) is similarly oriented for both triangles. The height of the subtriangle's vertex 1 above its opposite edge is h =û· (r p − r 1 ), p = 2 or 3 , and the (signed) coordinates of the two vertices at the FIGURE 2. Subtriangle T q with associated geometrical parameters. Note that T q can be further subdivided about the u-axis into two right sub-subtriangles. edge endpoints are (h, L , 0) and (h, U , 0), where L = ·(r 2 −r 1 ), U =ˆ ·(r 3 −r 1 ). Note that points along edges 2 and 3 (i.e., the edges opposite vertices 2 and 3 , resp.) where the principal branch of the arctangent is assumed. We note that, in the local subtriangle coordinate system, r − r = −uû− ˆ +dn and R = √ where d is a constant parameter and u and vary linearly over the subtriangle and parameterize source points on T q ; hence the static scalar and vector potential integrals (2) for general constant or linear bases may be constructed as superpositions of the following three subtriangle potentials with constant and linear source densities: for q = 1, 2, 3. Note that, for simplicity, we have excluded the usual material parameters , µ and 1/(4π ) factors in defining the integrals (3). Closed form expressions for the above integrals are summarized in the following three subsections. In each case, the subtriangle potential integral is integrated first with respect to (over the range u tan φ L ≤ ≤ u tan φ U ), and then with respect to u (on the range 0 ≤ u ≤ h). The result of the first, inner integral on (in square brackets following the second equality in (9)-(11) below) is given for completeness and verifiability. The result of the second, outer integration can be obtained by automatic symbolic integration, but usually requires manual simplification; all results have been checked by analytical differentiation to see that they reduce to the square bracketed integrands in (9)- (11). For compactness, the following quantities are defined and used in the expressions to follow (c.f., Fig. 2): (4) 99808 VOLUME 8, 2020 We also make use of the following defined auxiliary functions with removable singularities at x = 0: where sinhc −1 x follows the naming convention of other well-known cardinal functions, notably sinc x = sin x x and sinhc x = sinh x x . In addition, we also define Both sinhc −1 x and g(x) are even in x, non-negative, and infinitely differentiable, but to maintain accuracy, their approximate power series forms should be used to evaluate them near their removable singularity, i.e. for |x| < 10 −2 .
For both functions, all eight series terms should be retained for quadruple precision; for double precision applications, only the first five terms are needed. We also note that both functions are bounded for all real x and vanish as follows as |x| → ∞ : For convenience in plotting in the d = 0 constant plane, we also choose to evaluate the following indeterminate forms as shown:

A. ELECTROSTATIC POTENTIAL, UNIT SOURCE
The potential for subtriangle T q with a constant (unit) source density is given by where In arriving at (9), two arctangent terms appearing in I q (u, , d) are combined using the identity arctan . Forms equivalent to (9) have previously appeared in the literature [1], [16], but here the identity sinh −1 x = ln ( has been used to replace the logarithmic forms commonly arising from symbolic or tabular integration in favor of the more compact and convenient sinh −1 x function representation. Use of sinh −1 x, now an intrinsic function in most modern scientific programming languages, avoids loss in precision that can occur in the logarithmic form for large or small values of x. (The potential of a line source of length 2x observed at a distance y away in a plane passing through its center is sinh −1 (x/y); hence this form's frequent appearance in potential integral expressions should not be surprising.) In (9) and the following integrals, important limits as the observation point approaches either the u-plane of the subtriangle are also given. Testing for and correctly handling these limiting cases in (9) and the similar situations that follow below should be incorporated in any algorithm designed to evaluate them.

B. ELECTROSTATIC POTENTIAL, SOURCE LINEAR IN u
The potential for a subtriangle T q with a linearly varying source of density u is given by where Note use of the sinhc −1 x function in (10) to conveniently handle a removable singularity that occurs at P = 0, i.e., along lines normal to and passing through vertices of T q .

C. ELECTROSTATIC POTENTIAL, SOURCE LINEAR IN
The potential of subtriangle T q for a linearly varying source of density is given by where Here again, the sinhc −1 x function in (11) conveniently treats the removable singularity at P = 0.

D. STATIC VERTEX POTENTIAL FUNCTIONS
Consider the term I (u, , d) defined following Eq. (9). We may also write it as I (u, , d)| = u=u; =0 , in the notation of (9) since, at the lower limit, I (u, 0, d) = 0; thus I (u, , d) also represents the potential at (0, 0, d) of the unit source density, right sub-subtriangle shown in Fig. 3a with vertices at (0, 0, 0), (u, 0, 0), (u, , 0). Under this interpretation, the last term of Eq. (9) represents a superposition of potentials due to the two sub-subtriangles formed by splitting T q into two subsubtriangles about the u-axis. Similar observations apply to terms I u (u, , d) and I (u, , d), except that both have linear source densities, and, in the latter case, I (u, 0, d) = 0. 1 In any case, all three terms may be associated with potentials at the fixed observation point (0, 0, d) and having constant or linear source densities on the sub-subtriangle with vertices at (0, 0, 0), (u, 0, 0), and (u, , 0), relative to the observation point. Alternatively, if we merely translate the coordinate system origin to the latter vertex (leaving source distributions unchanged), as in Fig. 3b, the potentialĨ (x, y, z) in (x, y, z)coordinates at the same observation point in the new system is clearlyĨ (−u, − , d). But we also haveĨ (−u, − , d) = I (u, , d) = I (−u, − , d) since, for the first equality, the potential at any observation point is invariant with respect to coordinate translation, and, for the second, I (u, , d) is an even function of both u and . Similar arguments apply to I u (u, , d) and I (u, , d), except the latter are odd functions of u and , respectively. In summary, Thus, to within a sign factor and for fixed d, I , I u , and I may be interpreted either as (1) potentials at the observation point with u and locating the position of the source subsubtriangle vertex common to T and T q (Fig. 3a), or as (2) the potentials at the observation point whose projection is at −u and − with respect to the common vertex (Fig. 3b). In the first case, u and play the roles of locating source points with respect to the observation point; in the second, −u and − locate the observation point with respect to the source point. This simple consequence of the translational invariance of coordinate systems, that to within a sign factor, changes in u or for fixed d may be interpreted either as changes in source  or observation point locations, will be useful in determining potential gradients in Section 3.
We note that, in general, a linear combination of six terms of the form (9), (10), or (11) are needed altogether to define the (static) potential of a parent triangle for linear source distributions, and each is essentially independent of the parent triangle's shape. In general, each depends only on (a) the orientation of its associated parent triangle edge and (b) the position of its associated vertex along that edge, both located relative to the observation point.
We consider next an important property of singularities of the vertex potential functions defined following (9), (10), and (11). As partial static scalar potentials for sources defined on right triangles, vertex functions are analytic everywhere except in the source plane d = 0; there they must at least be continuous. Furtheremore, any low-order singularities along a parent triangle edge can arise solely from the pair of vertex potentials associated with that edge. And for observation points not directly on a parent edge, but rather on its extension beyond the parent's vertices, singularities from the vertex function pair must exactly cancel to maintain the analyticity of potentials along the source-free extended region. Since the vertex function pair are also independent, and the cancellation must also occur independent of the separation between an edge's vertices, the coefficient of any singularity must be constant along the extended edge. Thus an individual vertex's contribution is doubled between edges where they are additive, while it is cancelled along extended edge regions beyond the vertices. We note also that any singular behavior near a vertex arises solely from the two vertex functions defined there, each associated with a different local edge orientation.
A primary use of vertex functions is to reveal potential behavior near edges or vertices, independent of the shape of the source triangle T . We remark in passing that merely replacing the 1 R kernel in the vertex function integral definitions with the corresponding dynamic form, exp(−jkR) R extends the concept to time-harmonic potentials, with the the static case under current discussion corresponding to the special case k = 0, and providing the dominant contribution to the time-harmonic result. The potential of vertex functions for analyzing singular behavior for test integral quadrature design is briefly explored in the Appendix.

III. GRADIENTS OF SUBTRIANGLE POTENTIALS FOR CONSTANT AND LINEAR SOURCE DENSITIES
In this section, we evaluate the corresponding (vector-valued) integrals for gradient kernels in a somewhat unconventional way. We first note that if ψ is a potential quantity at an observation point due to a given source distribution on a parent triangle T (or sub-subtriangle), its value remains unchanged as long as the relative displacements between the observation point and vertices of the source (sub-)triangle are unchanged (i.e., the potential and source sub-subtriangle together are translationally invariant). Indeed, by the gradient property, a small relative displacement δr of the observation point with respect to the right source sub-subtriangle with vertices at (0, 0, 0), (u, 0, 0), (u, , 0) results in the potential change δψ ∼ = δr · ∇ψ(u, , d). In particular, since d is the relative distance alongn between the source triangle and observation point, a small displacement δd results in the potential change δψ ∼ = δdn · ∇ψ or, for infinitesimal changes, ∂ d ψ = n · ∇ψ, where ∂ d represents partial derivative with respect to d. On the other hand, the same change δψ results if, instead, we displace the vertices (i.e., the source points) of a sub-subtriangle by a small displacement δr = −δr in the opposite direction, i.e., δψ ∼ = −δr ·∇ψ. Thus, since u and fix the vertex locations of a right sub-subtriangle, changes in its potential contribution with respect to changes in u or are given by δψ ∼ = −δuû·∇ψ or δψ ∼ = −δ ˆ ·∇ψ, or for infinitesimal changes, byû·∇ψ orˆ ·∇ψ, respectively. Strictly speaking, we should also displace the sub-subtriangle's vertex at the projection point r 0 , but as already noted, the final potential (and hence also its gradient) must be independent of the projection point. Having determined components of the gradient along three orthogonal directions, we can now write the gradient operator for a given subtriangle potential as with gradient potential contribution The gradient operator (13) applies, of course, also to the potentials I q u , and I q l . We note in passing that (13) may also be derived by converting potentials to the alternate forms on the right hand side of (12), recognizing that in the converted forms, coordinates u, , and d all now represent observation points, and hence that the gradient form ∇ =û∂ u +ˆ ∂ +n∂ d , applies [17]. After then using (12) to convert back to the original potential forms we again arrive at (13) and (14). Applying the gradient operator (13) and (14) in turn to each of the subtriangle potential quantities (9)- (11), results in ∇I q , ∇I q u , and ∇I q as defined and discussed below.

A. GRADIENT OF ELECTROSTATIC POTENTIAL, UNIT SOURCE
The potential gradient contribution for a constant (unit) source density on subtriangle T q is found as follows: where In a software implementation of (15), the component quantitiesû · ∇I q ,ˆ · ∇I q , andn · ∇I q may be computed directly from (15) above, and these quantities conveniently absorb the negative signs appearing in the transverse components of the gradient operator. Also note that for d = 0, as r (hence also r 0 ) approaches an edge of the parent triangle, the ucomponent of the subtriangle's gradient becomes unbounded, whereas its -component remains bounded; more generally, in the plane of the parent triangle, derivatives of potential transverse to edges (i.e.,) are bounded, whereas derivatives normal to edges are logarithmically singular. Of particular note is that the contribution to the derivative alongn has a jump discontinuity at d = 0 arising from the sgn(d) factor. Furthermore, the normal derivative contribution of each subtriangle is proportional to the interior vertex angle at r(= r 0 ), tan −1 U /h − tan −1 L /h. When r is strictly interior to the parent triangle T , this angle sum over the three subtriangles is 2π , which, when the factors sgn(d) and 1/(4π ) are included, has the value ∓1/2 as d → 0 ± , a well-known factor (times the source density there, which here is unity) that usually appears outside the integral in integrals equations involving gradient kernels. As r approaches a triangle edge, however, one of the three subtriangles disappears completely, the angle sum becomes π , and after dividing by 4π , yields a coefficient of ∓1/4 on the source term instead. For r at a vertex of the parent triangle, only a single subtriangle contributes, however, and its (4π -normalized) contribution is just ∓α q /(4π ), VOLUME 8, 2020 where the parent (and subtriangle's) interior vertex angle is α q . Finally, the sgn(d) terms vanish since the angle sum is zero if r(= r 0 ) is strictly outside T . These properties of the potential gradient near a planar source are in agreement with the classical observations of [4], [18], [19] and others. Note that since the normal gradient potential is discontinuous across a source, merely specifying d = 0 in a procedure call is insufficient: Well-designed software should also provide a flag allowing one to specify which side of the surface source one wishes to observe, d = 0 + (top) or d = 0 − (bottom); merely specifying d = 0 should either result in an error message to the user and program termination or a warning message informing the user what default assumption has been made. Also note that as an observation point approaches a vertex of the parent triangle from an arbitrary direction, some gradient components approach limits that depend on the approach angles, θ and φ (c.f. Fig. 2). For example, in the plane d = 0, ∂ I q (u, , d) in (15) approaches u/P = cos φ, so that the limit at u = = 0 depends on the direction from which one approaches the point; lacking this information, the value of u/P there is indeterminate. We adopt the pragmatic view that well-designed quadrature schemes should avoid sampling the potential at such points, and hence we assign a value there, mostly for plotting purposes only; a logical choice for lim P→0 u/P, for example, would be its angular average over points in the infinitesimal disk in the plane d = 0 and centered at P = 0. Thus the line reading '' u/P P=0 := 0 '' in (15) implies that, although lim P→0 u/P is indeterminate at, P = 0 we assign it the value zero there. Finally, we observe in passing that removing the static part from a dynamic kernel can be especially important in computing potential gradients, as it essentially isolates all anomalous behavior to terms in (15), leaving the potential gradient integral with a difference kernel more suitable for numerical integration [15].

B. GRADIENT OF ELECTROSTATIC POTENTIAL, SOURCE LINEAR IN u
The potential gradients for sources linear in either u or are also computed by applying the gradient operator ∇ = −û∂ u −ˆ ∂ +n∂ d . For a source linear in u, we have the following result: where Note the use of the functions sinhc −1 (P/d) and g(P/d) to automatically handle removable singularities that would otherwise appear for P = 0 (for any d), i.e., for observation points on any of the three lines normal to the source plane and passing through vertices of the parent triangle T .

C. GRADIENT OF ELECTROSTATIC POTENTIAL, SOURCE LINEAR IN
The potential gradient for a source density linear in is as follows: where Again, the cardinal functions sinhc −1 (P/d) and g(P/d) automatically handle removable singularities at P = 0. We also observe that since the source densities u and both vanish at r 0 , the vertex gradient potentials ∇I q u of the previous subsection and ∇I q of this subsection are neither discontinuous nor unbounded for any r, including at edges or vertices.

IV. SYNTHESIS OF PARENT TRIANGLE POTENTIALS AND GRADIENTS FROM SUBTRIANGLE CONTRIBUTIONS
The total static potential of a parent triangle T with a unit source density is a superposition of its subtriangle potentials, (18) and the potential gradient is where I q and ∇I q are defined in (9) and (15), respectively. Combining (15) and (19), it is easily verified that the normal component of the gradient potential (19) as T is approached from above or below has the well-known classical limit [4], [18], [19]: where α(r) is the angular extent interior to T of a small disk centered about r = r 0 , i.e., 2π, r (= r 0 ) interior to T , π, r (= r 0 ) on an edge of T , α q , r (= r 0 ) at vertex q of T , 0, r (= r 0 ) exterior to T . (21) Note that the normal gradient component is discontinuous at d= 0, and hence in a software implementation, as previously noted, the intended limit for d= 0 can only be reliably guaranteed via a flag passed by the user. Similar superpositions of subtriangle potentials apply to potentials for linear source densities as well, i.e., (18) and (19) remain valid if we merely add subscripts u or to I , I q and their gradients. Beginning in this section and for the remainder of the paper, we restore the subtriangle index q as a superscript to subtriangle quantities that should be distinguished by their association with specific subtriangles. We observe, for instance, that a given subtriangle's potential contributions I q , I q u , or I q to the parent triangle potential may generally be neglected if |h q | < ε| q U − q L | where ε is a small (user-defined) positive number, say on the order of ε = 10 −14 for double precision results. This observation does not apply to the gradient potentials, however, since a vanishing potential does not imply the same of its gradient.
Turning next to the (static) vector potential integral, we note that for r in subtriangle T q , the pth RWG basis function (i.e., that associated with edge or vertex p of T ) can be expressed as for r in T ; it vanishes otherwise. Hence for the corresponding vector potential, we have and for its curl, via the identity ∇×(Aψ) = ∇ψ ×A+ψ∇× A, This expression satisfies the following classical limit as r approaches T : or equivalently, where α(r) is defined in (21) and the subscript ''tan'' indicates the transverse component.

V. NUMERICAL RESULTS
In this section numerical results are presented in the form of plots in the z = 0 plane of the static scalar potential and its gradient components for a unit density source on an equilateral triangle also lying in the z = 0 plane. For the same equilateral triangle and observation plane, plots of components of the static curl of vector potential for an RWG basis function source density are presented next. The section concludes with a series of plots illustrating the characteristic behaviors of the various vertex potentials and their gradients. As fundamental building blocks for all static potentials VOLUME 8, 2020 and their gradients for constant and linear source densities on planar elements, their discontinuities and singularities are important for understanding the behavior of the various potential quantities. Figure 4 (a) shows the scalar potential I (r) of (18) in the plane z = 0 of an equilateral triangle T with vertex locations (1, 0, 0) and (−1/2, ± √ 3/2, 0) and a unit density source. The vertex locations are easily identified in the figure, but edge locations, where slopes normal to an edge are infinite, are more difficult to visually resolve. They are better resolved in Fig. 4 (b), in which the vertical dimension of the plot is logarithmically compressed, but are best seen by plotting the tangential gradient componentsx·∇I andŷ·∇I as in Figs. 5 and 6, respectively, where they are singular at edges. Figure  5(a) plotsx·∇I while Fig. 5(b) shows a cut along the line y = −0.375, z = 0 and clearly displays the (logarithmically) singular nature of the gradient component normal to an edge as the edge is approached both from the interior and the exterior of T . By contrast, as Fig. 6 illustrates forŷ · ∇I ,   3/2 is the height of the equilateral triangle T . Note thatx · (∇ × I 1 ) is discontinuous at z = 0; values for z = 0 + are negatives of those for z = 0 − . gradient components parallel to an edge are bounded. As can be seen from the figures, if both highly accurate and efficient testing procedures for integral equations are to be employed, one must properly deal with both weak and logarithmically singular potentials and their gradients near edges. 99814 VOLUME 8, 2020 FIGURE 9. Plot of the y -componentŷ·(∇×I 1 ) =x· 1 (r)α(r)/(4π)

A. POTENTIALS AND GRADIENT POTENTIALS FOR A UNIT SOURCE EQUILATERAL TRIANGLE
3/2 is the height of the equilateral triangle T . Note thatŷ· (∇ × I 1 ) is discontinuous at z = 0. One possibilty is to use Duffy transforms to map test triangle domains to a rectangular domain, after which appropriate dyadic products of Gauss-Legendre and MRWlog integration rules can be used [20].  show potential integrals for which the integrands effectively involve the normal component of ∇(1/R). As is well-known and suggested by (20) and (21), such integrals are discontinuous as the observation point crosses the plane of the surface source, with discontinuity proportional to both the source density and the fractional portion of an infinitesimal disk at the observation point that lies in the interior of the source domain there. Figure 7 illustrates this for the normal component of the potential gradient, just on the underside (z = 0 − ) of same (unit) source triangular domain as in Figs. 4-6. On the topside (z = 0 + ), the normal gradient potential is merely the negative of that shown in Fig. 7. Figures 8 and 9 show tangential components of ∇×I p of (24) for p = 1, representing contributions to tangential magnetostatic field components just below the surface (z = 0 − ) of an RWG basis function, 1 (r), associated with vertex 1 (or equivalently, edge 1 opposite) on the x-axis of the equilateral triangle. The x-and y-components of the integral's curl are −∂ z I 1y and ∂ z I 1x , respectively, and in the plane d = 0 − the associated z-derivatives simply result in sampling the source function there, according to (21) and (22), and as reported in   Figs. 8 and 9. As Fig. 10 illustrates, the normal component of ∇ × I 1 is singular at domain boundaries.

B. VERTEX POTENTIAL FUNCTIONS AND THEIR GRADIENTS
In this section we explore plots of the vertex potentials and their gradients. Figure 11 shows a plot of I (u, , d) in the   plane d = 0. We note the function is bounded everywhere, but has an infinite slope along the line u = 0, d = 0, representing the (extended) edge of its associated source domain. As can be seen, the vertex potential I (u, , d) is odd in both variables u and ; it is also even in d. Furthermore, adding a subscript u or to I (u, , d), (i.e., I (u, , d) → I α (u, , d), α = u   or ) or differentiating I (u, , d), I u (u, , d), or I (u, , d) with respect u, , or d reverses the parity of the functions, as the plots following Fig. 11 verify. Taken together, these properties imply that it is sufficient to know a vertex potential or gradient potential function in a single octant, say u, , d > 0, in order to reconstruct it everywhere.

VI. CONCLUSION
Straightforward derivations of the classical potential functions and their gradients for constant and linear source distributions on triangles are presented. The various representations involve indeterminant and discontinuous forms as well as removable, bounded, and unbounded singularities. Hence, special consideration is given to representing results in computationally convenient form and to the proper evaluation or assignment of meaningful values to various limiting cases. Both Fortran and C++ coded versions of the resulting algorithms are available from the authors upon request. The concept of vertex potential functions and their gradients is introduced. Each function may be associated with a vertex and an (extended) edge of a planar, polygonal source domain. With the exception of possible jump discontinuities normal to the source plane, all anomalous behavior of a vertex function is confined to its associated vertex or extended edge. The potential for constant or linear source distributions on a triangle, for example, is a superposition of six vertex potentials. The domain-independence and singularity isolation properties of vertex potential functions are expected to prove useful in developing integration schemes for potentials and their gradients near triangle source boundaries occurring in the solution of integral equations by boundary element methods. This potential is briefly explored in the Appendix.
Numerical results are presented in graphical form illustrating the properties of potentials and potential gradients near

APPENDIX
To briefly illustrate the potential use of vertex functions in designing appropriate testing schemes for moment methods, consider the scalar potential reaction integral between a source and test triangle pair sharing a common edge. We assume the planes containing the source and test triangles make an angle β between their intersection along the -axis with origin at one of the common vertices, as shown in Fig. 21. Our focus is on modeling the potential behavior over the test triangle in the neighborhood of that vertex. We select the -axis as the polar axis of an (R, α, β) spherical coordinate system centered at the vertex, where R is the radial distance from the common vertex to the observation point and α is its angle with respect to the -axis. In terms of spherical coordinates, the (u, , d) coordinates of the observation point with respect to the source triangle vertex are thus given by u = R sin α cos β = R 0 cos β, = R cos α, d = −R sin α sin β, where R 0 = √ u 2 + d 2 = √ R 2 − 2 = R sin α and tan α = R 0 / . Thus, near the vertex, the dominant potential vertex function contribution to the scalar potential (9), in spherical coordinates, becomes I q (u, , d) = R sin α cos β sinh −1 (cot α) − |sin β| tan −1 cos α cos β sin α + |sin β| .
We see immediately that, for fixed β, the R and α dependencies are conveniently separated, and the test triangle integral may therefore be written as cos α cos β sin α + |sin β| dα, VOLUME 8, 2020 where α U is the vertex angle of the test triangle at its common vertex with the source triangle, and R(α) expresses the angular dependence of the upper limit for radial integration. A more detailed definition for R(α) depends either on the test triangle shape or (more likely) how one decides to split the test triangle into subtriangles to confine test integrals to neighborhoods about each common vertex. In any case, for the remaining angular integral, we see that the term sinh −1 (cot α) is logarithmically singular along the polar axis where α = 0 or π. The factor sin α also vanishes at α = 0 and π and hence reduces the order of the integrand singularity; nevertheless, special handling will still be required to achieve high accuracy in the testing integral. Similar considerations apply to the other potential integrals or their gradients.