Generalized Debye Sources Based EFIE Solver on Subdivision Surfaces

TL;DR

This paper introduces a novel, well-conditioned GDS-based EFIE solver on subdivision surfaces, combining generalized Debye sources with isogeometric analysis for improved accuracy and efficiency in electromagnetic scattering problems.

Contribution

It extends the GDS method to subdivision surface discretizations, enabling well-conditioned systems and efficient isogeometric analysis for electromagnetic scattering.

Findings

Demonstrates improved conditioning of the EFIE system.

Shows high accuracy and efficiency of the GDS-IGA framework.

Validates the approach with numerous numerical results.

Abstract

The electric field integral equation is a well known workhorse for obtaining fields scattered by a perfect electric conducting (PEC) object. As a result, the nuances and challenges of solving this equation have been examined for a while. Two recent papers motivate the effort presented in this paper. Unlike traditional work that uses equivalent currents defined on surfaces, recent research proposes a technique that results in well conditioned systems by employing generalized Debye sources (GDS) as unknowns. In a complementary effort, some of us developed a method that exploits the same representation for both the geometry (subdivision surface representations) and functions defined on the geometry, also known as isogeometric analysis (IGA). The challenge in generalizing GDS method to a discretized geometry is the complexity of the intermediate operators. However, thanks to our earlier work on subdivision surfaces, the additional smoothness of geometric representation permits discretizing these intermediate operations. In this paper, we employ both ideas to present a well conditioned GDS-EFIE. Here, the intermediate surface Laplacian is well discretized by using subdivision basis. Likewise, using subdivision basis to represent the sources, results in an efficient and accurate IGA framework. Numerous results are presented to demonstrate the efficacy of the approach.

Figures

Equations

$$\hat{n}\times \mathcal{T}(\mathbf{j}({\mathbf{r}}^{\prime }))=-\hat{n}\times \hat{n}\times {\mathbf{E}}^{inc}(\mathbf{r}),\mathbf{r}\in \Gamma ,$$

$$\mathcal{T}(\mathbf{X}(\mathbf{r}))=\hat{n}\times ik\eta {\int }_{\Gamma }\left[\mathcal{I}+\frac{\nabla \nabla }{k^2}\right]g(\mathbf{r},{\mathbf{r}}^{\prime })\cdot \mathbf{X}({\mathbf{r}}^{\prime })d{\mathbf{r}}^{\prime }$$

$$g(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{e^{ikR}}{4\pi R}$$

$$\begin{array}{rl}\mathbf{j}(\mathbf{r})& ={\nabla }_{\Gamma }\Psi (\mathbf{r})+\hat{n}\times {\nabla }_{\Gamma }\Phi (\mathbf{r})+{\mathbf{j}}_H(\mathbf{r})\\ & ={\mathbf{j}}_1(\mathbf{r})+{\mathbf{j}}_2(\mathbf{r})+{\mathbf{j}}_H(\mathbf{r}),\end{array}$$

$${\nabla }_{\Gamma }\cdot {\mathbf{j}}_H(\mathbf{r})=0,{\nabla }_{\Gamma }\cdot \left({\hat{n}\times \mathbf{j}}_H(\mathbf{r})\right)=0.$$

$${\nabla }_{\Gamma }^2\Psi (\mathbf{r})={\Delta }_{\Gamma }\Psi (\mathbf{r})=i\omega p(\mathbf{r}),$$

$${\nabla }_{\Gamma }^2\Phi (\mathbf{r})={\Delta }_{\Gamma }\Phi (\mathbf{r})=-i\omega q(\mathbf{r}).$$

$${\mathbf{j}}_1(\mathbf{r})=i\omega {\nabla }_{\Gamma }{\Delta }_{\Gamma }^{-1}p(\mathbf{r}),{\mathbf{j}}_2(\mathbf{r})=-i\omega \hat{n}\times {\nabla }_{\Gamma }{\Delta }_{\Gamma }^{-1}q(\mathbf{r}),$$

$$\hat{n}\times \mathcal{T}\left(\mathcal{C}(p({\mathbf{r}}^{\prime }))-{\hat{n}}^{\prime }\times \mathcal{C}(q({\mathbf{r}}^{\prime }))\right)=-\hat{n}\times \hat{n}\times {\mathbf{E}}^{inc}(\mathbf{r}),\mathbf{r}\in \Gamma ,$$

$$\mathcal{C}(f(\mathbf{r}))=i\omega {\nabla }_{\Gamma }{\Delta }_{\Gamma }^{-1}f(\mathbf{r}).$$

$${\mathcal{M}}_{\Gamma }=\left\{f:\Gamma \to \mathbb{R}|{\int }_{\Gamma }f(\mathbf{r})d\mathbf{r}=0\right\}.$$

$${\nabla }_{\Gamma }^2\Psi (\mathbf{r})=f(\mathbf{r}),$$

$${\int }_{\Gamma }\Psi (\mathbf{r})d\mathbf{r}=0.$$

$$\mathbf{S}(\mathbf{r}(u,v))=\sum _{n=1}^{12}{\mathbf{t}}_n{\xi }_n(\mathbf{r}(u,v))$$

$$p(\mathbf{r})\approx \sum _{n=1}^N{a}_{1,n}{\xi }_n(\mathbf{r}),q(\mathbf{r})\approx \sum _{n=1}^N{a}_{2,n}{\xi }_n(\mathbf{r}),$$

$$F(\Psi )=\frac{1}{2}{\int }_{\Gamma }{\nabla }_{\Gamma }^2\Psi (\mathbf{r})d\mathbf{r}-{\int }_{\Gamma }f(\mathbf{r})\Psi (\mathbf{r})d\mathbf{r}$$

$$\mathbf{G}{\mathbf{b}}_1=\mathbf{h}$$

$$G_{m,n}={\int }_{{\Gamma }_m\cap {\Gamma }_n}{\nabla }_s{\xi }_m(\mathbf{r})\cdot {\nabla }_s{\xi }_n(\mathbf{r})d\mathbf{r}$$

$$h_n=-{\int }_{{\Gamma }_n}{\xi }_n(\mathbf{r})f(\mathbf{r})d\mathbf{r}.$$

$$\sum _{i=1}^N{b}_{1,n}{c}_n=0.$$

$$\mathbf{T}=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 0& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1\\ -c_1/c_N& -c_2/c_N& \cdots & -c_{N-1}/c_N\end{array}\right]$$

$$\mathbf{G}\mathbf{T}{\mathbf{b}}_1^{\prime }=\mathbf{h}.$$

$${\mathbf{T}}^T\mathbf{G}\mathbf{T}{\mathbf{b}}_1^{\prime }={\mathbf{T}}^T\mathbf{h}.$$

$$P(\Psi )=F(\Psi )+\frac{\gamma }{2}{\left({\int }_{\Gamma }\Psi (\mathbf{r})d\mathbf{r}\right)}^2$$

$${\mathbf{G}}^p{\mathbf{b}}_1=\mathbf{h}$$

$$\begin{array}{rl}{G}_{m,n}^p=& {\int }_{{\Gamma }_m\cap {\Gamma }_n}{\nabla }_s{\xi }_m(\mathbf{r})\cdot {\nabla }_s{\xi }_n(\mathbf{r})d\mathbf{r}\\ & +\gamma {\int }_{{\Gamma }_m}{\xi }_m(\mathbf{r})d\mathbf{r}{\int }_{{\Gamma }_n}{\xi }_n(\mathbf{r})d\mathbf{r}.\end{array}$$

$$L(\Psi ,\beta )=F(\Psi )+\beta {\int }_{\Gamma }\Psi (\mathbf{r})d\mathbf{r}.$$

$$\left[\begin{array}{cc}\mathbf{G}& \mathbf{c}\\ {\mathbf{c}}^T& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{b}}_1\\ \beta \end{array}\right]=\left[\begin{array}{c}\mathbf{h}\\ 0\end{array}\right].$$

$$\Psi (\mathbf{r})\approx \sum _{n=1}^N{b}_{1,n}{\xi }_n(\mathbf{r}),\Phi (\mathbf{r})\approx \sum _{n=1}^N{b}_{2,n}{\xi }_n(\mathbf{r}).$$

$${\mathbf{j}}_1(\mathbf{r})\approx \sum _{n=1}^N{a}_{1,n}{\mathbf{j}}_{1,n}(\mathbf{r})=i\omega \sum _{n=1}^N{b}_{1,n}{\nabla }_{\Gamma }{\xi }_n(\mathbf{r}),$$

$${\mathbf{j}}_2(\mathbf{r})\approx \sum _{n=1}^N{a}_{2,n}{\mathbf{j}}_{2,n}(\mathbf{r})=-i\omega \sum _{n=1}^N{b}_{2,n}\hat{n}\times {\nabla }_{\Gamma }{\xi }_n(\mathbf{r}).$$

$$\left[\begin{array}{cc}{\mathbf{Z}}^{11}& {\mathbf{Z}}^{12}\\ {\mathbf{Z}}^{21}& {\mathbf{Z}}^{22}\end{array}\right]\left[\begin{array}{c}{\mathbf{b}}_1\\ {\mathbf{b}}_2\end{array}\right]=\left[\begin{array}{c}{\mathbf{e}}_1\\ {\mathbf{e}}_2\end{array}\right]$$

Full text

I Introduction

Surface integral equation (SIE) solvers have been the mainstay in computational electromagnetics, including a range of problems [1, 2, 3]. In particular, the electric field integral equation (EFIE) has been extensively explored for various applications [4, 5]. As a result, considerable research effort has been invested in understanding the nuances and ramifications of discretizing these equations. This includes various efforts to understand low frequency breakdown [6, 7, 8], develop well conditioned formulations [9], introduce higher order basis sets [10] and investigate accuracy and convergence [11], develop hierarchical basis [12], and so on. However, by and large, the problem has remained the same: how does one develop integral formulations that are well behaved across frequencies of interest especially when high discretization density is required to capture geometric features. To this end, several new SIE formulations and numerical techniques have been proposed; a partial listing of these includes the current and charge integral equation (CCIE) [9], augmented EFIE (A-EFIE) [13], Calderón preconditioner [14], multi-resolution analysis [15] and introducing loop-tree/star basis functions [7, 8, 16, 17, 18], and Debye sources [19, 20, 21]. Recently, a decoupled potential integral equation (DPIE) based on Lorentz gauge has been proposed in [22] that leads to a second-kind and stable formulation over a wide frequency band. DPIE was implemented numerically by using Nyström method in [23]. A similar idea of exploiting generalized gauge based ${\bf A}$-$\Phi$ integral formulation was introduced in [24] and its numerical implementation was presented in [25].

The focus of this work is to build on an approach presented by Epstein and Greengard [19]. Their approach relies on using two scalar sources and building an SIE framework to solve for these unknowns. In their work, they demonstrated the efficacy of the approach when applied to a sphere and later to arbitrary shapes. The crux of this approach is to define two scalars, called Debye sources that are employed to represent Debye potentials which, in turn, can be generalized to represent currents on the arbitrary surface, not limited to the sphere surface. Two second-kind scalar integral equations can be derived by using this framework. However, implementing these equations in a discrete setting is challenging as one needs to find the inverse of the surface Laplacian or the Laplace-Beltrami operator that maps the unknown Debye sources to Debye potentials. A recent paper by Chernokozhin and Boag [26] presented the first numerical discrete implementation of [19] on piecewise smooth surfaces. In their work, the inverse of surface Laplacian operator was obtained by exploiting the finite difference scheme defined on structured grids on the piecewise flat surface. Besides, additional constraints have to be imposed in order to ensure the continuity of surface currents. The question we seek to ask is whether one can create a better/more accurate framework on surface representation that are described by an underlying local parameterization.

To summarize the state of art of Debye sources based integral equation solvers, methods have been presented for analytical surface descriptions, and piecewise constant/lowest order Lagrangian surface description. The principal challenge is the lack of a robust surface Laplacian equation solver. The approach that we present overcomes this bottleneck by using a smooth definition of the surfaces. To this end, we take recourse to our recent work on isogeometric methods [27]. Iso-geometric/parametric analysis (IGA) has seen a resurgence in recent years, thanks in large part to Hughes et al. [28]. The rationale for this research has been to provide a seamless interface between computer aided design (CAD) and computer aided engineering (CAE). In this framework, the same basis function set is employed to represent both geometry and underlying physical quantities residing on the geometry. While a bulk of the existing work has focused on using variations of splines, in what follows, we use subdivision surface representations [29] that provide $C^{2}$ smoothness on the geometry almost everywhere. This representation was exploited to develop an IGA based EFIE solver [27] that relied on a surface Helmholtz decomposition, wherein the current is represented in terms of surface gradient and surface curl of potentials; henceforth referred to as P-EFIE. In this paper, we will show that this additional smoothness provides the means to obtain a convergent inverse of the surface Laplacian without imposing additional constraints to ensure continuity.

The principal focus of this work is to exploit isogeometric basis sets to extend the idea of generalized Debye sources (GDS) [19, 20, 21] to traditional EFIE. This formulation will be referred to as GDS-EFIE. We will employ subdivision basis sets to create a discrete system to solve surface Laplacian equation embedded in GDS-EFIE first and then resolve the whole GDS-EFIE system within the IGA framework. In this paper, we shall demonstrate

• •

performance of different numerical techniques for solving surface Laplacian equation,

• •

convergence of the eigenvalues estimation and the solution to the inverse of Laplace-Beltrami operator,

• •

low frequency stability and better conditioning of the IGA GDS-EFIE solver,

• •

flexibility of the solution technique.

Parenthetically, we note that the proposed solver is straightforward to implement compared to the formulations proposed in [19, 26].

The remainder of this paper is organized as follows: Section II presents the classical EFIE briefly and then proposes the variation of the EFIE formulation based on generalized Debye sources. Subdivision surfaces and functions are presented concisely in Section III. Numerical implementations of GDS-EFIE including surface Laplacian equation solvers are detailed in Section IV. Numerical examples validating the proposed approach are demonstrated in Section V. Finally, Section VI summarizes the contribution of this work as well as future directions.

II Formulations

II-A Electric Field Integral Equation

Consider a simply-connected PEC object $\Omega$ whose boundary $\partial\Omega$ is denoted by $\Gamma$ which is equipped with a unique outward pointing normal $\hat{n}$. This object is illuminated by incident fields (${{\bf{E}}^{inc}}$, ${{\bf{H}}^{inc}}$), and it results in scattered fields exterior to $\Omega$. Total fields denoted by $({\bf E},{\bf H})$ are composed of incident fields and scattered fields. The scattered fields due to the object can be obtained by using the equivalence principle, which posits an equivalent current ${\bf j}({\bf r})$ for ${\bf r}\in\Gamma$ such that the total magnetic field just outside the surface $\Gamma^{+}$ satisfies the jump condition ${\hat{n}}\times{\bf H}={\bf j}$, with a slight abuse of notation. The total electric field satisfies ${\hat{n}}\times{\bf E}=0$, which can be rewritten as

[TABLE]

that formulates the electric field integral equation (EFIE). In the above equation, the integral operator ${\cal{T}}$ is defined as

[TABLE]

in which $i$ is the imaginary unit ($i\equiv\sqrt{-1}$), ${k}$ is the wavenumber, ${\eta}$ is the intrinsic impedance of the background medium, ${\cal{I}}$ is the identity operator and

[TABLE]

is the Green’s function in free-space. In the above equation, ${R=|{\bf{r}}-{\bf{r}}^{\prime}|}$ is the distance between the field point ${\bf{r}}$ and the source point ${\bf{r}}^{\prime}$. An $e^{-i\omega t}$ time-dependence convention is assumed and suppressed throughout this paper. Here, $\omega$ is the angular frequency. The solution to the EFIE is typically effected by using method of moments (MoM) wherein surface current is represented by a set of vector basis functions, say the Rao-Wilton-Glisson basis functions [30] which are equivalent to the lowest order Raviart-–Thomas functions [31]. Alternatives to this approach has been a topic of significant recent interest; these include using generalized method of moments (GMM) [32, 33], subdivision surfaces [27], discontinuous basis set [34] , and more recently, Debye sources [19]. All the aforementioned methods try to bring features into modeling electromagnetic scattering; but a common thread that ties GMM, MoM on subdivision surfaces, and Debye sources is the use of surface Helmholtz decomposition. In what follows, we prescribe that generalized Debye sources can be exploited to solve the EFIE. Ideas for using Debye sources were initiated by Epstein and others in a series of papers [20, 19, 21]. While our approach is slightly different, genesis of ideas are rooted in the above citations. This approach necessitates the use of scalar representations for the sources. Next, we briefly describe EFIE based on generalized Debye sources.

II-B Generalized Debye Sources based EFIE

Consider an arbitrary current ${\bf j}(\bf r)$ that exists on surface $\Gamma$. It is well known that this current may be written by using a surface Helmholtz decomposition as

[TABLE]

where $\nabla_{\Gamma}$ is the surface gradient, $\Psi({\bf r})$ and $\Phi({\bf r})$ are scalar potentials and ${\bf j}_{H}({\bf r})$ is a harmonic vector field that satisfies

[TABLE]

The harmonic component vanishes (${\bf j}_{H}({\bf r})\equiv 0$) on the surface of simply-connected geometry. The potentials can be related to two scalar sources $p({\bf r})$ and $q({\bf r})$, the so-called generalized Debye sources via:

[TABLE]

It is noted that $p({\bf r})$ has direct physical meaning here, i.e., surface charge density unlike the one in [19, 26], and $\Delta_{\Gamma}$ is the surface Laplacian or Laplace-Beltrami operator. Then the currents can be rewritten as

[TABLE]

in which $\Delta^{-1}_{\Gamma}$ is the inverse of Laplace-Beltrami operator. As surface current is represented by generalized Debye sources, EFIE (1) can be rewritten in terms of $p({\bf r}^{\prime})$ and $q({\bf r}^{\prime})$ as follow:

[TABLE]

where

[TABLE]

Equation (8) can be considered as generalized Debye sources based EFIE (GDS-EFIE). It is noted that there are always non-trivial solutions to equation (6), $\Psi$=constant and $\Phi$=constant. To ensure uniqueness of the inverse of Laplace-Beltrami operator, $\mathcal{M}_{\Gamma}$, a space of mean-zero function on $\Gamma$ can be defined first:

[TABLE]

Then surface Laplacian operator $\Delta_{\Gamma}$ might be invertible if a map from $\mathcal{M}_{\Gamma}$ to itself on $\Gamma$ can be established [19]. Given that function space of generalized Debye sources belongs to $\mathcal{M}_{\Gamma}$, a unique solution to the inverse of $\Delta_{\Gamma}$ can be achieved by solving the following boundary value problem:

[TABLE]

Equation (11b) can be considered as a constraint to surface Laplacian equation (11a) in order to obtain a unique solution. Similar comments are valid for $\Phi({\bf r})$, and are omitted in the above discussion for the sake of brevity.

It is worth noting that both discretization and seeking the inverse of the surface Laplacian operator are non-trivial. First, one needs to find a smooth enough representation for each scalar quantity on complex surfaces so that the resulting currents satisfy the continuity condition. The work in [26] introduced additional conditions to guarantee continuity. Our proposed method, however, relies on a smooth basis set ($C^{2}$ almost everywhere and $C^{1}$ at isolated points) for scalars with the help of the subdivision surface representation. Second, the direct inverse of surface Laplacian operator doesn’t exist if no constraint is imposed to remove the known constant-value null space. In the next section, subdivision surface and subdivision basis set are reviewed briefly.

III Subdivision Surfaces and Functions

As a shape description, subdivision surfaces technique has been explored extensively in computer graphics, especially the animation industry. Even though non-uniform rational B-splines (NURBS) is built in most CAD systems as an industry standard, NURBS results in a smooth description in the interior of a patch but only $C^{0}$ or even worse across the boundary between patches [35]. In comparison, the limit surface generated by subdivision schemes is $C^{2}$ almost everywhere except at finite points (irregular vertices) where the description is $C^{1}$. Thus, one can exploit subdivision basis to resolve GDS-EFIE within IGA framework without imposing any additional conditions to ensure continuity.

Since triangular tessellations are omnipresent in SIE solvers for EM problems, we employ function spaces defined by the Loop subdivision scheme [29]. Consider an initial primal mesh denoted by $P_{0}$ at level 0 which consists of a set of vertices and connectivity map. The valence of a given vertex is defined as the number of triangles incident on itself. A 1-ring of a vertex consists of all vertices of these triangles. A vertex is considered as a regular one if its valence equals to 6. Otherwise, it is called an irregular or extraordinary vertex. A triangle is regular if its vertices are all regular, and irregular otherwise. A limit surface can be generated by recursive refinements of primal mesh $P_{0}$. Specifically, there are $\mathcal{N}_{t}$ triangles in $P_{0}$ and Loop subdivision scheme generates triangular tessellations of level $k$ recursively including $4^{k}\mathcal{N}_{t}$ new triangular patches by inserting new vertices at the edge midpoints (in parameter domain sense) and subdividing a triangle of level $k-1$ into four sub-triangles of level $k$. It is noted that every newly inserted vertices are regular and only irregular vertices at level 0 still remain irregular. After each subdivision, the position of every vertex will be recomputed and each new triangle patch can be parameterized by 1-ring of the patch(union of vertices of incident triangles to its three vertices). Fig. 1 illustrates a regular triangle $\mathcal{E}$ defined by its 1-ring vertices indexed from 1 to 12.

When a nodal value function $\xi_{n}({\bf r})$ is associated with $n$th vertex, a regular triangular patch can be evaluated as

[TABLE]

where $(u,v)$ denotes a pairwise coordinate on a parameterization chart and ${\bf t}_{n}\in\mathbb{R}^{3}$ is the position vector of $n$th vertex. For a regular triangular patch, $\xi_{n}({\bf r})$ is a quartic box spline tabulated in [36]. Since a triangular patch is defined by its 1-ring vertices, the scalar function $\xi_{n}({\bf r})$ can affect on 2-ring domain of vertex ${\bf t}_{n}$ (union of the 1-rings of the 1-ring). It is zero outside the vertex’s 2-ring. Thus, $\xi_{n}({\bf r})$ has a compact support and $\xi_{n}({\bf r})\in C^{2}_{0}$. For an irregular vertex, $\xi_{n}({\bf r})$ is still a piecewise polynomial function, however, it has a infinite number of segments towards the irregular vertex itself. As such, $\xi_{n}({\bf r})$ will degrade to $C^{1}$ smoothness only at the irregular vertex. For more details on evaluating $\xi_{n}({\bf r})$ associated with an irregular vertex, we refer to the reference [36, 27].

IV Algorithmic Flowchart

In this section, we will describe the procedure that we will exploit to solve GDS-EFIE (8). The solution to this equation proceeds from right to left. Specifically, the unknowns in the system are the scalar sources, $p({\bf r})$ and $q({\bf r})$. These are then mapped to potentials, $\Psi({\bf r})$ and $\Phi({\bf r})$, which are in turn mapped to the equivalent electric current ${\bf j}({\bf r})$. Finally, the integral operator ${\cal T}$ maps the current onto the scattered electric field. In what follows, each of these mapping operations is discussed in sequence, starting with the inverse of the surface Laplacian and followed by representation of current and then the discretization of the operator.

To start, we assume that given $\xi_{n}({\bf r})$ for $n$th vertex of the primal mesh $M_{0}$ with $N$ vertices in total, the unknown generalized Debye sources $p({\bf r})$ and $q({\bf r})$ in GDS-EFIE can be represented by

[TABLE]

where $a_{1,n}$ and $a_{2,n}$ are unknown coefficients. Since surface current representations in terms of generalized Debye sources require the inverse of Laplace-Beltrami operator, three different numerical techniques for solving surface Laplacian equation are discussed next. In what follows, we only illustrate the solution for $\Psi({\bf r})$; the solution for $\Phi({\bf r})$ can be obtained in a similar manner.

IV-A Surface Laplacian Solvers

It can be shown that the solution of equation (11a) will minimize the functional

[TABLE]

which can be achieved by imposing the stationary requirement ${\delta}F=0$. To satisfy this requirement, we start with representing the potential $\Psi({\bf r})$ by using subdivision basis functions such that $\Psi({\bf r})\approx\sum\nolimits_{n=1}^{N}b_{1,n}\xi_{n}({\bf r})$. Then a $N\times N$ matrix equation corresponding to the weak form of variational problem $\delta F=0$ can be written as:

[TABLE]

where

[TABLE]

Here, $\Gamma_{m(n)}$ denotes the support of basis function $\xi_{m(n)}$. Thus far, only discretization of equation (11a) is considered. Since matrix ${\bf G}$ is rank deficient by one due to the existing of one-dimensional null space, one needs to solve the equation (11a) in tandem with (11b). In this work, three different techniques including least squares (LSQ), penalty and Lagrange multiplier (LM) are studied. Comparisons in effectiveness and efficiency are also made among these methods in the results section.

IV-A1 Least Squares Solution

Substituting the representation of the potential function to the constraint (11b) and letting $\int_{\Gamma_{n}}{\xi_{n}}({\bf r})d{\bf r}=c_{n}$, we have

[TABLE]

Therefore, the last coefficient $b_{1,N}$ can be expressed in a linear combination of others as $b_{1,N}=-{{\bf b}^{\prime}_{1}}^{T}{\bf c}^{\prime}/c_{N}$, where ${\bf b}^{\prime}_{1}=(b_{1,1},b_{1,2},...,b_{1,N-1})^{T}$ and ${\bf c}^{\prime}=(c_{1},c_{2},...,c_{N-1})^{T}$. As a result, one can define a sparse transformation matrix

[TABLE]

that maps ${\bf b}^{\prime}_{1}$ onto ${\bf b}_{1}$ via ${\bf T}{\bf b}^{\prime}_{1}={\bf b}_{1}$. The constraint (11b) together with the discrete system can be rewritten as

[TABLE]

This is an overdetermined $N\times(N-1)$ matrix equation with respect to ${\bf b}^{\prime}_{1}$ and may be solved approximately via least squares method as follow:

[TABLE]

It is noted that the product of ${\bf T}^{T}{\bf G}{\bf T}$ is not necessarily sparse. However, as each involved matrix of the product is sparse, one may obtain an iterative solution with ${\cal O}(N)$ for per matrix vector multiplication (MVM).

IV-A2 Penalty Method

While the above approach takes the constraint (11b) directly into account, an alternate approach is augmenting the functional with a penalty constraint. The functional defined earlier may be modified to include a penalty function as follows:

[TABLE]

where $\gamma$ is a prescribed parameter. Hence, the matrix equation (15) might be modified as

[TABLE]

where

[TABLE]

Here, $p$, the subscript of $G_{m,n}^{p}$, indicates the matrix entry of penalty method. Similar to LSQ method, the final system is a full matrix.

IV-A3 Lagrange Multiplier Method

Another viable alternative to the above is to employ a Lagrange multiplier; as before, this involves a minor modification to the functional that is to be minimized, viz.,

[TABLE]

Similarly, the stationary requirement $\delta L=0$ results in a unique ${\bf b}_{1}$ and $\beta_{0}$ that will minimize the above functional. Letting $\partial L/\partial b_{1,n}=0$ and $\partial L/\partial\beta=0$, one can obtain a matrix equation with an augmented dimension of $(N+1)\times(N+1)$:

[TABLE]

It is noted that the resulting matrix of Lagrange multiplier method is sparse when $N\gg 1$. Hence, the complexity for both memory and MVM scales as $\mathcal{O}(N)$.

IV-B Discretization of GDS-EFIE

Thus far, we have prescribed three different approaches for seeking the inverse of the surface Laplacian operator based on the following representation:

[TABLE]

It follows that once the inverse is found, we can obtain the coefficients $b_{1,n}$ and $b_{2,n}$ which are related to $a_{1,n}$ and $a_{2,n}$ that are involved to represent the unknown Debye sources $p({\bf r})$ and $q({\bf r})$. Then, surface currents ${\bf j}_{1}({\bf r})$ and ${\bf j}_{2}({\bf r})$ can be expressed as

[TABLE]

Using these expressions for the current and testing GDS-EFIE (8) with $\nabla_{\Gamma}\xi_{n}({\bf r})$ and $\hat{n}\times\nabla_{\Gamma}\xi_{n}({\bf r})$ results in a system that reads as [27]

[TABLE]

where

[TABLE]

Then an iterative solution for the unknown coefficients $a_{i,n}$ for $i=1,2$ proceeds via two sets of iterations (i) inverse of the surface Laplacian to map from $a_{i,n}\longrightarrow b_{i,n}$ for $i=1,2$, and (ii) convergence of (28). Together, this results in solutions to $a_{i,n}$. Next, we discuss a series of results that address convergence of surface Laplacian solvers as well as solutions to the proposed EFIE.

V Numerical examples

V-A *Performance Comparison of Surface Laplacian Solvers *

To test the performance of three aforementioned surface Laplacian solvers, three metrics are employed: (a) the residual of left hand side of constraint (11b) which is supposed to be zero exactly, (b) the condition number of the resulting matrix system and (c) the convergence rate for each solver. For a given right hand side of (11) $f\in\mathcal{M}_{\Gamma}$, matrix equations (20), (22) and (25) are solved iteratively with maximum of iteration count set to 600 and tolerance $10^{-13}$, respectively. Bi-conjugate gradients (BiCG) method is employed as an iterative solver for surface Laplacian equation. In this paper, condition number of matrix ${\bf A}$ is defined in matrix 2-norm sense, namely,

[TABLE]

For penalty method, the prescribed parameter $\gamma$ ranges from 1 to $10^{9}$ and is sampled at every order. Finally, the examples chosen are akin to those used in computer graphics [37].

A sphere with radius $r=1$m is first analyzed. The sphere surface is discretized into 5120 triangular elements. Hence, there are 2562 vertices. The residual of constraint of penalty method as a function of parameter $\gamma$ is depicted in Fig. 2. For better comparison and illustration, the residuals of least square (LSQ) and Lagrange multiplier (LM) methods are also plotted in the figure even though they are both independent of $\gamma$. As observed, the solution of penalty method is more accurate as $\gamma$ increases. It is also evident that the accuracy of both LSQ and LM method is very high.

The condition number of the matrix system of three different solvers is plotted in Fig. 3, respectively. From the figure, it is evident that the condition number of penalty method increases rapidly as $\gamma$ increases. The condition number of LM method is the lowest among three solvers.

Fig. 4 depicts the relative residuals of three solvers as iteration number increases. For penalty method, two curves corresponding to $\gamma=1$ and $\gamma=10^{7}$ respectively are plotted.

In the second example, we consider a hand shape object meshed as shown in inset of Fig. 5. The hand is discretized by 2880 triangular elements and 1442 vertices. Fig. 5 demonstrates the residuals of constraint (11b) and it is evident that the accuracy of penalty method is the worst. Even when $\gamma=10^{9}$, the residual of penalty method is two orders larger than LSQ and LM method. Again, as parameter $\gamma$ increases, the condition number of penalty method will be larger as shown in Fig. 6. LM method outperforms once more. Finally, Fig. 7 depicts the residual changes as iteration number increases for three solvers. It is noted that the accuracy of penalty method is poor when $\gamma=1$ although it can achieve the tolerance fast. Hence, it is not straightforward to provide a guideline how to choose parameter $\gamma$ for penalty method. The accuracy and matrix conditioning will go to the opposite side when $\gamma$ increases.

Taking the three metrics into account, we posit that the performance of LM method is optimal amongst the three solvers. Moreover, the implementation of LM is not expensive since the complexity of per MVM only scales as $\mathcal{O}(N)$.

V-B Convergence Test of the Surface Laplacian

While the above results presented how well different methods worked, it is also important that the results by themselves converge with spatial refinement. To this end, we formulate a generalized eigenvalue problem and exploit this as a metric to study convergence, specifically that of the LM method. The generalized eigenvalue can be formulated as

[TABLE]

where $\lambda_{n}$ is $n$th generalized eigenvalue to solve, ${\bf v}$ is the corresponding eigenvector and

[TABLE]

Here, ${\bf M}$ is a mass matrix with its entry defined as

[TABLE]

We solve the generalized eigenvalue problems for the both sphere and hand shape object analyzed earlier. To study the experimental order of convergence, an initial control mesh with 642 vertices for the sphere is used and three subsequent subdivisions are conducted which can generate 2562, 10242 and 40962 vertices, respectively. The reference solution comes from one more subdivision based on the finest mesh and there are 163842 vertices in total. Fig. 8 plots the relative errors in three generalized eigenvalues $\lambda_{400}$, $\lambda_{500}$ and $\lambda_{600}$ for a sphere as a function of mesh size $h$. It indicates the convergence rate scales as $h^{2p}$ approximately where $p=3$ when subdivision basis function is employed. As for hand shape object, we start from an initial control mesh with 1442 vertices. Three subsequent subdivisions are processed as well corresponding to 5762, 23042 and 92162 vertices respectively. Again, the reference solution is computed based on one more subdivision with 368642 vertices. The relative error of three generalized eigenvalues $\lambda_{1200}$, $\lambda_{1300}$ and $\lambda_{1400}$ is demonstrated in Fig. 9. From these results, it is apparent that the surface Laplacian equation system can be well discretized by using subdivision basis $\xi({\bf r})$ within IGA framework.

V-C EM Scattering from PEC Objects

Next, several numerical examples are presented to illustrate the efficacy of the proposed GDS-EFIE for analysis of EM scattering problems. In all cases, we present comparison with of radar cross section (RCS) data obtained using the proposed method with that against a well validated Rao-Wilton-Glisson based EFIE code, and against analytical data (when available).

First, we consider a PEC sphere with radius $r=0.67\lambda$ illuminated by a plane wave traveling along $-\hat{z}$ with electric field $-\hat{x}$ polarized is considered first. Here $\lambda$ is the working wavelength. There are 642 vertices for discretization. Bistatic RCS solutions from GDS- and P-EFIE [27] are plotted in Fig. 10. Both EFIE formulations have an excellent agreement with the reference solution obtained by Mie series method. Fig. 11 depicts the relative residuals when the iterative solver is applied to both GDS- and P-EFIE. With generalized Debye sources as unknowns, GDS-EFIE converges faster compared with P-EFIE.

Next, to demonstrate the low-frequency stability of GDS-EFIE, the mesh size is fixed and the working frequency decreases gradually. Fig. 12 plots the iteration numbers to achieve the prescribing tolerance $10^{-5}$ for GDS-EFIE for frequencies ranging from 1Hz to 10MHz. In the figure, the frequency is sampled at every order. We can see that the iteration number is stable in the low-frequency regime. It is noted that P-EFIE will not converge without any preconditioning techniques.

The second example involves scattering from a warhead which fits into a box $10.4\lambda\times 3.5\lambda\times 3.5\lambda$. The incident plane wave propagates along $-\hat{z}$ with electric field $-\hat{x}$ polarized. The warhead is discretized by using 21376 triangle elements and then there are 10690 vertices. Fig. 13 depicts the surface current density of the warhead obtained by GDS-EFIE and it is evident that there are no artificial defects since smooth subdivision basis function is applied. Bistatic solutions from subdivision basis for GDS-EFIE and RWG based EFIE code are plotted in Fig. 14. It is evident that both solutions agree with each other very well.

As the last example, scattering from a plane model with electrical size of $6.6\lambda\times 6.6\lambda\times 1.8\lambda$ is analyzed. The object is illuminated by a plane wave traveling along $\hat{y}$ with electric field $\hat{z}$ polarized. The plane is meshed by 39984 triangle elements associated with 19994 vertices. Fig. 15 shows the surface current density and again, there does not exist any artifacts. Bistatic solutions are plotted in Fig. 16 and solution of GDS-EFIE agrees well with the reference solution from RWG based EFIE.

VI Conclusion

In this work, we have developed an integral formulation called GDS-EFIE based on scalar unknowns by extending the idea of generalized Debye sources to traditional EFIE. The proposed formulation inherits the salient property of the work in [19, 26], i.e., well-conditioned integral equation at low frequency regime. As is evident in the above description, GDS-EFIE is straightforward to implement numerically. The challenge lies on solving the surface Laplacian equation both effectively and efficiently. Thanks to the smooth subdivision basis, we can find the inverse of Laplace-Beltrami operator within IGA framework with high order accuracy and convergence. Furthermore, no additional continuity constraint is required and our framework can be applied to arbitrarily shaped simply-connected triangular control mesh, without being limited to flat structured mesh. Several numerical examples have been presented to show the well conditioning and flexibility of our proposed GDS-EFIE at both regular and low frequency regimes. Our next steps are to develop a solver for the original Debye source method proposed in [19] by using subdivision basis sets within IGA framework, and then extend these to composite/multiply connected objects.

Acknowledgment

This work was supported in part by Hong Kong GRF 712612E, NSFC 61271158, US AOARD 124082 and 134140 contracted through UTAR, and Hong Kong UGC AoE/P04/08. The authors thank the high performance computing center at Michigan State University for computational support.