An isogeometric boundary element method for electromagnetic scattering with compatible B-spline discretizations

TL;DR

This paper presents an isogeometric boundary element method using compatible B-splines for electromagnetic scattering, enabling high-order accuracy and direct CAD integration for complex geometries.

Contribution

It introduces a novel compatible B-spline construction for electromagnetic boundary element analysis that directly utilizes CAD data, improving efficiency and accuracy.

Findings

Successfully verified with Mie scattering and NASA almond problems.

Achieves efficient computation with H-matrices and Bézier extraction.

Handles complex geometries directly from CAD without meshing.

Abstract

We outline the construction of compatible B-splines on 3D surfaces that satisfy the continuity requirements for electromagnetic scattering analysis with the boundary element method (method of moments). Our approach makes use of Non-Uniform Rational B-splines to represent model geometry and compatible B-splines to approximate the surface current, and adopts the isogeometric concept in which the basis for analysis is taken directly from CAD (geometry) data. The approach allows for high-order approximations and crucially provides a direct link with CAD data structures that allows for efficient design workflows. After outlining the construction of div- and curl-conforming B-splines defined over 3D surfaces we describe their use with the electric and magnetic field integral equations using a Galerkin formulation. We use B\'ezier extraction to accelerate the computation of NURBS and B-spline terms and employ H-matrices to provide accelerated computations and memory reduction for the dense matrices that result from the boundary integral discretization. The method is verified using the well known Mie scattering problem posed over a perfectly electrically conducting sphere and the classic NASA almond problem. Finally, we demonstrate the ability of the approach to handle models with complex geometry directly from CAD without mesh generation.

Figures

Tables

Table 1. Table 1 : Details of div-conforming B-spline discretizations used in the Mie scattering study.

mesh (# elements)	degrees of freedom
	$(1,0)\times (0,1)$	$(2,1)\times (1,2)$	$(3,2)\times (2,3)$	$(4,3)\times (3,4)$
h0 (6)	12	48	108	192
h1 (24)	48	108	192	300
h2 (96)	192	300	432	588
h3 (384)	768	972	1,200	1,452
h4 (1536)	3,072	3,468	3,888	4,332

Table 2. Table 2 : Details of compatible B-spline discretizations used for the NASA almond study.

mesh (# elements)	degrees of freedom
	$(1,0)\times (0,1)$	$(2,1)\times (1,2)$	$(3,2)\times (2,3)$
h0 (288)	558	700	858
h1 (1152)	2,268	2,546	2,840
h2 (4608)	9,144	9,694	10,260

Equations

$$\mathbf{n}\times \mathbf{E}=0$$

$$\mathbf{E}={\mathbf{E}}^i+{\mathbf{E}}^s$$

$${\mathbf{E}}^s=-j\omega \mathbf{A}-\nabla \phi$$

$$\phi (\mathbf{x})=\frac{1}{\epsilon }{\int }_{\Gamma }\rho \frac{e^{-jkr}}{4\pi r}\mathrm{d}\Gamma (\mathbf{y})$$

$$\rho =-\frac{1}{j\omega }\nabla \cdot \mathbf{J}$$

$$\mathbf{A}(\mathbf{x})=\mu {\int }_{\Gamma }\mathbf{J}(\mathbf{y})\frac{e^{-jkr}}{4\pi r}\mathrm{d}\Gamma (\mathbf{y}).$$

$${\mathbf{E}}^s=-j\omega \mu \left({\int }_{{\Gamma }_y}\mathbf{J}\frac{e^{-jkr}}{4\pi r}\mathrm{d}{\Gamma }_y+\frac{1}{k^2}{\nabla }_{{\Gamma }_x}{\int }_{{\Gamma }_y}{\nabla }_{{\Gamma }_y}\cdot \mathbf{J}\frac{e^{-jkr}}{4\pi r}\mathrm{d}{\Gamma }_y\right)$$

$$L^E[\mathbf{\tau }(\mathbf{x})]={\int }_{{\Gamma }_y}\mathbf{\tau }\frac{e^{-jkr}}{4\pi r}\mathrm{d}{\Gamma }_y+\frac{1}{k^2}{\nabla }_{{\Gamma }_x}{\int }_{{\Gamma }_y}{\nabla }_{{\Gamma }_y}\cdot \mathbf{\tau }\frac{e^{-jkr}}{4\pi r}\mathrm{d}{\Gamma }_y$$

$$⟨\mathbf{w},L^E[\mathbf{J}]⟩=⟨\mathbf{w},\mathbf{f}⟩\forall \mathbf{w}\in \mathcal{V}$$

$$⟨{\mathbf{w}}_h,L^E[{\mathbf{J}}_h]⟩=⟨{\mathbf{w}}_h,\mathbf{f}⟩\forall {\mathbf{w}}_h\in {\mathcal{V}}_h.$$

$$\mathbf{F}=\sum _{a=1}^{n_p}{R}_a(\mathbf{s}){\mathbf{P}}_a$$

$$R_a(\mathbf{s})\equiv R_a(s,t)=\frac{w_a{B}_a(s,t)}{{\sum }_{b=1}^{nm}{w}_b{B}_b(s,t)}a=1,2,\dots n_p$$

$$B_a(s,t)=B_i^p(s)B_j^p(t),$$

$$a=(j-1)n+i.$$

$${\mathcal{M}}_h=\{Q=({\zeta }_i^1,{\zeta }_i^1+1)\times ({\zeta }_j^2,{\zeta }_j^2+1),1\le i\le n_k-1,1\le j\le m_k-1\}$$

$${\hat{S}}^p:=\text{span}\{B_i^p{\}}_{i=1}^n$$

$${\hat{S}}^{p,q}:={\hat{S}}^p\otimes {\hat{S}}^q=\text{span}\{B_a{\}}_{a=1}^{n_b}$$

$${\hat{S}}_1:={\hat{S}}^{p,q-1}\times {\hat{S}}^{p-1,q}$$

$${\hat{S}}_2:={\hat{S}}^{p-1,q}\times {\hat{S}}^{p,q-1}.$$

$${\mathcal{U}}_h=\{\mathbf{u}:\mathbf{u}\circ \mathbf{F}=\frac{1}{J}D\mathbf{F}\hat{\mathbf{v}},\hat{\mathbf{v}}\in {\hat{S}}_1\}$$

$${\mathcal{V}}_h=\{\mathbf{v}:\mathbf{v}\circ \mathbf{F}={\left(D{\mathbf{F}}^{+}\right)}^{\mathrm{T}}\hat{\mathbf{v}},\hat{\mathbf{v}}\in {\hat{S}}_2\}$$

$$D\mathbf{F}=\left[\begin{array}{cc}\frac{\partial x}{\partial s}& \frac{\partial x}{\partial t}\\ \frac{\partial y}{\partial s}& \frac{\partial y}{\partial t}\\ \frac{\partial z}{\partial s}& \frac{\partial z}{\partial t}\end{array}\right],$$

$$D{\mathbf{F}}^{+}={\left(D{\mathbf{F}}^{T}D\mathbf{F}\right)}^{-1}D{\mathbf{F}}^T,$$

$$J=\sqrt{{\left(\frac{\partial y}{\partial s}\frac{\partial z}{\partial t}-\frac{\partial z}{\partial s}\frac{\partial y}{\partial t}\right)}^2+{\left(\frac{\partial z}{\partial s}\frac{\partial x}{\partial t}-\frac{\partial x}{\partial s}\frac{\partial z}{\partial t}\right)}^2+{\left(\frac{\partial x}{\partial s}\frac{\partial y}{\partial t}-\frac{\partial y}{\partial s}\frac{\partial x}{\partial t}\right)}^2}.$$

$${\hat{\mathbf{N}}}_a^{\text{div}}(s,t)=\{\begin{array}{l}{B}_a^{(p,q-1)}(s,t){\hat{\mathbf{e}}}_{1}1\le a\le n_b^1\\ B_{a-n_b^1}^{(p-1,q)}(s,t){\hat{\mathbf{e}}}_2{n}_b^1+1\le a\le n_b^1+n_b^2\end{array}$$

$${\mathbf{N}}_a^{\text{div}}(\mathbf{x}(s,t))=\frac{1}{J}D\mathbf{F}{\hat{\mathbf{N}}}_a^{\text{div}}(s,t)1\le a\le n_b=n_b^1+n_b^2$$

$${\mathbf{u}}_h^{\text{div}}(\mathbf{x})=\sum _{a=1}^{n_b}{\mathbf{N}}_a^{\text{div}}(\mathbf{x})u_a$$

$${\mathbf{v}}_h^{\text{curl}}(\mathbf{x})=\sum _{a=1}^{n_b}{\mathbf{N}}_a^{\text{curl}}(\mathbf{x})v_a$$

$$\overline{\Gamma }=\bigcup _{i=1}^{n_d}\overline{{\Gamma }_i}$$

$${\mathbf{F}}_i=\sum _{a=1}^{n_p^i}{R}_a^i(\mathbf{s}){\mathbf{P}}_a^i$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Research into unifying geometry and analysis for efficient design workflows has progressed rapidly in recent years driven by the isogeometric analysis and computational geometry research communities. Analysis based on geometry discretizations now covers a wide range of technologies including NURBS [1], T-splines [2], LR B-splines [3], PHT-splines [4] and subdivision surfaces [5]. A major research challenge at present is the automatic generation of volumetric discretizations from given geometric surface data and promising research includes the work of [6, 7] based on T-splines. In contrast, analysis methods based on shell formulations or boundary integral methods are known to require only a surface discretization exhibiting key benefits for a common geometry and analysis model since no additional volumetric processing is required. There has been much research into isogeometric shell formulations including [5, 8, 9] and developments into isogeometric boundary element methods based on NURBS [10, 11], T-splines [12, 13] and subdivision surfaces [14].

A key application of the boundary element method is the analysis of electromagnetic scattering over complex geometries in which a perfectly electrically conducting (PEC) assumption can be made. The method is often termed the method of moments within the electromagnetic research community but is synonymous with the Galerkin boundary element method. It is well known that a straightforward application of nodal basis functions to the electric and magnetic field integral equations (EFIE, MFIE) prevents numerical convergence and instead, discrete spaces that satisfy the relevant continuity requirements must be used. The most commonly used discretization that satisifes the relevant continuity requirements are Raviart-Thomas [15] or RWG [16] basis functions that are mainly based on low order polynomials.

In the context of isogeometric analysis progress has been made on the development of spline-based compatible discretizations [17, 18, 19, 20] in which a discrete de Rham sequence can be constructed providing a crucial step towards application of isogeometric analysis for fluid flow and electromagnetics applications. This fundamental work opens up the opportunity for the development of an isogeometric boundary element method (isogeometric method of moments) for electromagnetic scattering which is the focus of the present study. We note similar work in which subdivision surfaces are employed [21], but we believe that use of B-spline based algorithms provides greater refinement flexibility, provide a natural link with NURBS based systems that are ubiquitous in modern engineering design software, and offer higher convergence rates over equivalent subdivision schemes with extraordinary points.

We organise the paper as follows: first, we prescribe the Galerkin formulation of the relevant integral equations that govern electromagnetic scattering; we give an overview of NURBS surfaces and detail the construction of compatible B-splines; we then specify the fully discretized form of the integral equations for electromagnetic scattering with compatible B-splines; we cover implementation details of the method including fast evaluation of basis functions through Bézier extraction and the use of $\mathscr{H}$-matrices to approximate dense matrices; we verify the present method by performing electromagnetic scattering over a sphere in which a closed-form solution is provided by Mie scattering theory and finally, we demonstrate the ability of the present approach to perform electromagnetic scattering of PEC bodies with complex geometries taken directly from CAD software. It is assumed that time-harmonic fields are prescribed and, unless stated otherwise, it can be assumed that $\mathbf{x}\in\mathbb{R}^{3}$.

2 Electric field integral equation: Galerkin formulation

We first assume a PEC domain $\Omega$ with connected boundary $\Gamma:=\partial\Omega$ residing within an unbounded domain $\Omega_{\infty}$ with isotropic permeability and permittivity given by the scalar quantities $\varepsilon$ and $\mu$ respectively. We further assume a polarised time-harmonic electromagnetic plane wave of angular frequency $\omega$ is imposed on the PEC body with a wavenumber $k=\omega\sqrt{\varepsilon\mu}$. Denoting $\mathbf{E}$ as the total electric field, in the presence of an electromagnetic wave a surface current $\mathbf{J}$ is induced and the following PEC condition holds on the surface of the scattered object

[TABLE]

where $\mathbf{n}$ represents the outward pointing normal vector. We specify the incident wave as $\mathbf{E}^{i}(\mathbf{x})=\mathbf{p}\,e^{-jk\mathbf{d}\cdot\mathbf{x}}$ where $j$ is the unit imaginary number, $\mathbf{p}=(p_{x},p_{y},p_{z})$ is a polarization vector and $\mathbf{d}=(d_{x},d_{y},d_{z}),|\mathbf{d}|=1$ is a propagation vector. The relationship between the total, incident and scattered electric fields is written as

[TABLE]

where $\mathbf{E}^{s}$ represents the scattered electric field. The entire set-up is depicted in Figure 1.

Following the potential formulation of Maxwell’s equations (see e.g. [22]), the scattered electric field can be expressed in terms of an electric potential $\varphi$ and magnetic vector potential $\mathbf{A}$ (assuming time-harmonic fields) as

[TABLE]

where the electric potential is given by

[TABLE]

with $r:=|\mathbf{x}-\mathbf{y}|$ and the charge density $\rho$ expressed as

[TABLE]

with the magnetic potential related to the surface current through

[TABLE]

We omit variable dependencies in future equations where they are implied by their context and adopt the notation $\Gamma_{y}\equiv\Gamma(\mathbf{y})$ and $\Gamma_{x}\equiv\Gamma(\mathbf{x})$. Substituting (4) and (6) into (3) and employing (5) with $k^{2}=\omega^{2}\varepsilon\mu$ and $j^{2}=-1$, the scattered electric field is expressed in terms of surface quantites as

[TABLE]

where $\nabla_{\Gamma_{x}}$, $\nabla_{\Gamma_{y}}$ are surface gradient operators taken with respect to $\mathbf{x}$ and $\mathbf{y}$ respectively. Defining the linear operator

[TABLE]

along with the force term $\mathbf{f}=(j\omega\mu)^{-1}\mathbf{E}^{i}$, the Galerkin formulation of the EFIE reads as:

given $\mathbf{f}$, find $\mathbf{J}\in\mathcal{V}$ such that

[TABLE]

where $\mathcal{V}$ is the trace space $H^{-\frac{1}{2}}(\mathrm{div}_{\Gamma},\Gamma)$, and the $\langle\cdot,\cdot\rangle$ is the duality pairing between $\mathcal{V}$ and $H^{-\frac{1}{2}}(\mathrm{curl}_{\Gamma},\Gamma)$. When the fields are smooth enough, the duality pairing reduces to $\langle\mathbf{u},\mathbf{v}\rangle=\int_{\Gamma}\mathbf{u}\cdot\mathbf{v}\,\mathrm{d}\Gamma$.

We define the finite dimensional subspace $\mathcal{V}_{h}\subset\mathcal{V}$ which allows the solution of (9) to be approximated as the solution of

given $\mathbf{f}$, find $\mathbf{J}_{h}\in\mathcal{V}_{h}$ such that

[TABLE]

Conventionally, $\mathbf{w}_{h}$ and $\mathbf{J}_{h}$ are discretized through the Raviart-Thomas basis, but in our approach we make use of compatible B-splines that we now outline in detail.

3 Discretization

3.1 NURBS surfaces

Our implementation assumes a watertight NURBS surface parameterization that may be composed of multiple patches and we further assume that the connectivity of global basis functions between NURBS patches is known a priori. Dealing with the single patch case first, a NURBS surface parameterization is defined through a set of four-dimensional homogeneous control points $\{\mathbf{P}_{a}\}_{a=1}^{n_{p}}$, $\mathbf{P}_{a}=(x_{a}w_{a},y_{a}w_{a},z_{a}w_{a},w_{a})$ (where $w_{a}$ represents a control point weight), a set of knot vectors $\{\Xi_{i}\}_{i=1}^{2}$ where $\Xi_{1}=\{0=s_{1},s_{2},\ldots,s_{n+p+1}=1\}$, $\Xi_{2}=\{0=t_{1},t_{2},\ldots,t_{m+q+1}=1\}$ and a degree vector $\mathbf{p}=(p,q)$. $n$ and $m$ denote the number of basis functions defined through the knot vectors $\Xi_{1}$ and $\Xi_{2}$ respectively with $n_{p}=n\times m$. We assume all knot vectors are open (i.e. for a given degree $p$ the knot vector contains $p+1$ equal knot values at its beginning and end).

Defining the parametric domain $\widehat{\Gamma}=(0,1)^{2}\subset\mathbb{R}^{2}$ and physical domain $\Gamma\subset\mathbb{R}^{3}$, a NURBS geometric mapping $\mathbf{F}:\widehat{\Gamma}\to\Gamma$ can be written in terms of parametric coordinates $\mathbf{s}=(s,t)\in\widehat{\Gamma}$ as

[TABLE]

with the set of rational basis functions $\{R_{a}\}_{a=1}^{n_{p}}$ defined as

[TABLE]

where

[TABLE]

with the set of univariate B-spline basis functions $\{B^{p}_{i}\}_{i=1}^{n}$ defined through the Cox-de-Boor algorithm (see e.g. [23]). The parametric basis function index $a$ is defined in terms of the univariate basis indices $i,j$ through

[TABLE]

Defining vectors of unique knot values in the $s$ and $t$ parametric directions as $\boldsymbol{\zeta}_{1}=\{\zeta_{1}^{1},\zeta_{2}^{1},\ldots\zeta_{n_{k}}^{1}\}$ and $\boldsymbol{\zeta}_{2}=\{\zeta_{1}^{2},\zeta_{2}^{2},\ldots\zeta_{m_{k}}^{2}\}$ respectively, the mesh in the parametric domain is given by

[TABLE]

with $n_{e}=\textrm{size}(\mathcal{M}_{h})$ denoting the number of elements within the patch. Each element $Q$ within the patch contains $(p+1)\times(q+1)$ non-zero basis functions.

3.2 Compatible B-spline approximation

Given a set of univariate B-spline basis functions $\{B_{i}^{p}\}_{i=1}^{n}$, the space spanned by this basis is defined as

[TABLE]

and in a similar manner, the tensor product B-spline space defined through the set of B-spline basis functions $B_{a}:=B_{i}^{p}\otimes B_{j}^{q}$, $i=1,2,\ldots,n$, $j=1,2,\ldots m$ is defined as

[TABLE]

where the mapping defined by (13) is employed and a hat symbol denotes that the quantity is defined over the parametric domain. A div-conforming vector B-spline space is defined over the parametric domain as

[TABLE]

and likewise, a curl-conforming vector B-spline space is defined as

[TABLE]

The equivalent div-conforming and curl-conforming spaces defined in the physical domain are then constructed through appropriate Piola mappings as

[TABLE]

and

[TABLE]

respectively, where $D\mathbf{F}$ is the Jacobian associated with the geometric mapping $\mathbf{F}$ which for 3D surfaces is given by the rectangular matrix

[TABLE]

$D\mathbf{F}^{+}$ is the Monroe-Penrose pseudoinverse of the Jacobian given by

[TABLE]

and $J$ is the surface element given by

[TABLE]

Further details of the derivation of (19) and (20) can be found in [18, 20] and the derivation of (21)-(23) can be found in [24, Sect. 5.4].

3.2.1 Basis functions

Expressing vectors within the parametric domain as $\widehat{\mathbf{v}}=\widehat{v}_{i}\widehat{\mathbf{e}}_{i}$, $i=1,2$ and adopting the notation $\{B_{a}^{(p,q-1)}\}_{a=1}^{n_{b}^{1}}$, $\{B_{a}^{(p-1,q)}\}_{a=1}^{n_{b}^{2}}$ to represent the set of B-spline basis functions associated with the spaces $\widehat{S}^{p,q-1}$ and $\widehat{S}^{p-1,q}$ respectively, the set of div-conforming basis functions in the parametric domain $\hat{\Gamma}$ is defined as

[TABLE]

which are transformed into a set of div-conforming basis functions on the surface $\Gamma$ using the Piola transformation defined in (19) as

[TABLE]

where $\mathbf{F}\equiv\mathbf{F}(s,t)$ is implied. Curl-conforming basis functions are defined in analogous fashion.

Global div- and curl-conforming approximations in physical space can then simply be expressed through

[TABLE]

and

[TABLE]

respectively, where $u_{a}$ and $v_{a}$ are control coefficients. To illustrate the construction of compatible B-splines based on the NURBS parameterization shown in Figure 2, the bivariate B-splines generated from univariate B-splines are shown for two example basis functions in Figure 3. Further application of the Piola transformation as defined in (25) generates the div-conforming B-spline basis functions in physical space as shown in Figure 4.

Remark 1

For simplicity the construction of compatible B-splines is described using the same degree $(p,q)$ of the geometry. In practice it is possible to use a different degree for the B-splines discretization, as we will see in the numerical experiments.

3.3 Multipatch discretizations

Invariably, NURBS surfaces will consist of multiple patches whose union defines the physical domain through

[TABLE]

where $n_{d}$ is the number of parametric domains or patches and $\Gamma_{i}\cap\Gamma_{j}=\emptyset$ for $i\neq j$. Each domain $\Gamma_{i}$ is constructed through a NURBS geometric mapping $\mathbf{F}_{i}:\widehat{\Gamma}\to\Gamma_{i}$ with parametric coordinates $\mathbf{s}\in\widehat{\Gamma}$ as

[TABLE]

where the index $i$ indicates that the relevant quantity is restricted to patch $\Gamma_{i}$. We require for two patches $\Gamma_{i}$ and $\Gamma_{j}$ with $i\neq j$ and which share a common edge the geometry mapping along the shared edge is the same. In addition, the knot vectors associated with each patch at the common edge must be the same, up to an affine transformation. Figure 5(a) illustrates the geometry mappings of a multipatch NURBS geometry.

A global geometry connectivity array $C_{g}$ can be defined which maps a parametric basis function index $a$ and patch index $i$ to a global geometry basis index as

[TABLE]

The definition of the geometry connectivity array and the NURBS parameterisation given by (29) allows a multipatch NURBS parameterisation to be constructed such as that shown in Figure 5(b).

As is well-known with vector bases, care must be taken when constructing global compatible basis functions since both the global basis function index and the orientation sign must be stored and we refer the reader to [25] where div- and curl-conforming B-spline approximations are constructed in a volumetric context. We define the vector basis connectivity for a div-conforming basis through

[TABLE]

where $n_{b}^{i}$ is the number of compatible B-spline basis functions in patch $i$. This allows a global multipatch compatible B-spline discretization to be written as

[TABLE]

where $N_{b}$ is the global number of basis functions, $\mathbf{N}_{A}^{\textrm{div}}|_{\Gamma_{i}}\equiv\mathbf{N}^{\textrm{div}}_{C_{n}(i,a)}\equiv\textrm{sgn}(i,a)\mathbf{N}^{\textrm{div}}_{i,a}$.

From an implementation standpoint the main consideration is how to handle basis functions along the edges of parametric domains which is best illustrated graphically. Figure 6 shows an example vector basis connectivity for div-conforming B-splines of order $(4,3){\mkern-2.0mu\times\mkern-2.0mu}(3,4)$ based on the geometry of Figure 5. Similar connectivities can be constructed for curl-conforming B-splines.

4 Discretised EFIE with compatible B-splines

In the present work $\mathbf{w}_{h}$ and $\mathbf{J}_{h}$ in (10) are defined through the the div-conforming B-spline discretization given by (31) and can be expressed as

[TABLE]

Substituting (32) and (33) into (10) and applying the divergence theorem to transfer a derivative onto $\mathbf{w}_{h}$, a system of equations is formed as

[TABLE]

where

[TABLE]

[TABLE]

and $\mathbf{J}_{B}$ represents a vector of unknown surface current density coefficients. A similar procedure can be carried out for the magnetic field integral equation as detailed in Appendix A.

4.1 Radar Cross Section

The radar cross section $\sigma$ which quantifies how detectable an object is to a radar signal in a given direction is computed as

[TABLE]

where $R$ is the distance between the radar signal and the target object and furthermore, it can be assumed in the present work that $|\mathbf{E}^{i}|$ = 1. As detailed in [26, 27] if the source and field points are located far apart then $R\approx|\mathbf{x}|$ and the scattered electric field at a source (observation) point can be expressed as

[TABLE]

allowing the RCS to be computed as

[TABLE]

or, in terms of the RCS in decibels per square metre

[TABLE]

5 Implementation

Figure 7 details the main steps in the implementation of the present method. A multipatch compatible B-spline discretization is constructed directly from the NURBS surface parameterization. The inherent link between the geometry and analysis models allows for straightforward computation of compatible basis functions with the relevant Piola transforms. We utilise Bézier extraction [28] to accelerate computations whereby high order B-spline and NURBS basis functions are computed through precomputed Bézier extraction coefficients and inexpensive Bernstein polynomials.

As is well-known with Galerkin boundary element methods, careful consideration must be given to the computation of the matrix components $\mathbf{Z}_{AB}$ given by (35) when the element domains $\Gamma_{x}$ and $\Gamma_{y}$ are either coincident, edge adjacent, vertex adjacent or lie close to one another. We use the robust quadrature algorithms proposed by Sauter and Schwab [29] that deal with each of these cases.

To overcome the debilitating nature of large dense matrix $\mathbf{Z}$, we approximate this matrix using $\mathscr{H}$-matrices whereby a low-rank approximation is constructed through appropriate geometrical cluster trees that separate terms into admissible and non-admissible terms (i.e. far-field and near-field terms respectively). We do not wish to delve into the technical details of $\mathscr{H}$-matrices and instead guide the reader to relevant literature (see e.g. [30, 31]). However, we remark that $\mathscr{H}$-matrices are found to be particularly amenable for implementation into an existing BEM library and we make use of the library HLibPro [32] which provides high-performance $\mathscr{H}$-matrix libraries that scale optimally over multicore hardware and are primarily based on the Adaptive Cross Approximation algorithm [33]. The library requires as an input the set of bounding boxes defined by the support of each basis function (see Figure 8) and the basis function index associated with each box. Once an $\mathscr{H}$-matrix approximation is formed for a particular wavenumber, the matrix can be written and read freely from file which allows for highly efficient radar cross section computations. We note that this approach is valid for low to medium wavenumbers with special techniques required for high wavenumbers (e.g. [34]).

6 Numerical results

To verify the present approach and to demonstrate the capability of the method of performing electromagnetic scattering directly from CAD models using an isogeometric approach we present numerical results for a series of electromagnetic scattering problems with PEC conditions.

6.1 PEC sphere

The first problem we consider is that of electromagnetic plane wave impinging on a PEC sphere of radius $a=1$ which has a well-known solution given by the Mie series (see e.g. [22]). The incident wave is polarised in the x-direction by specifying $\mathbf{p}=(1,0,0)$ and is chosen to propagate in the positive z-direction with $\mathbf{d}=(0,0,1)$. The solution for the surface current given in spherical coordinates $(\rho,\theta,\phi)$ (see Figure 10) is expressed as

[TABLE]

with

[TABLE]

where $\eta=\sqrt{\mu/\varepsilon}$, the terms $P_{n}^{1}$ and $P_{n}^{1^{\prime}}$ correspond to the set of order 1 associated Legendre polynomials and derivatives respectively and

[TABLE]

with $h^{(2)}_{n}$ denoting the spherical Hankel function of the second kind. The radar cross section for this problem given in terms of increasing normalised wavenumber is illustrated in Figure 9 where the two asymptotic limits associated with Rayleigh and optical scattering are labelled.

Using the present approach, the sphere geometry is discretised using bi-quartic NURBS patches arranged in a cube topology with no degenerate points, as in Figure 5(a). Control point coordinates, weights and knot vectors for this NURBS parameterization can be found in [35]. We construct div-conforming B-splines using the knots inherited by the NURBS parameterization with degrees $(4,3){{\mkern-2.0mu\times\mkern-2.0mu}}(3,4)$, $(3,2){{\mkern-2.0mu\times\mkern-2.0mu}}(2,3)$, $(2,1){{\mkern-2.0mu\times\mkern-2.0mu}}(1,2)$ and $(1,0){{\mkern-2.0mu\times\mkern-2.0mu}}(0,1)$ and apply successive h-refinement (knot insertion) to generate a set of meshes h0 (base mesh), h1, h2 etc. Table 1 provides further details of each discretization. It should be noted that compatible B-splines of degree $(1,0){{\mkern-2.0mu\times\mkern-2.0mu}}(0,1)$ are directly equivalent to low order Raviart-Thomas or RWG basis functions on quadrilateral meshes. The bi-quartic NURBS representation of the geometry is used for all analyses and thus geometric error is eliminated for all discretizations considered.

After solving for surface current, equations (38) and (39) were used to determine radar cross section values with the results for mesh h3 shown in Figure 11 for each B-spline degree. The superior RCS accuracy obtained through higher order B-spline discretizations is demonstrated and this is also apparent in RCS values obtained with meshes h0, h1 and h2 as presented in Appendix B. As expected, finer meshes are capable of handling higher wavenumbers.

Plots of surface currents and magnitudes for $k=8$, h3 are shown for each B-spline degree in Figures 12 through to 15 where the higher accuracy and smoothness offered through higher B-spline degrees is visible. Recalling that the $(1,0){{\mkern-2.0mu\times\mkern-2.0mu}}(0,1)$ discretization is equivalent to the commonly used Raviat-Thomas elements, it is clear that higher order compatible B-spline discretizations offer substantial accuracy improvements over such elements.

Additionally, to establish that correct convergence rates are obtained using our approach we compute relative errors using the norm defined by

[TABLE]

where we remark that the $L^{2}$ norm of the surface divergence is well defined for this particular example. A convergence rate of $p+1$ is expected for a given B-spline degree with minimum degree $p$. We specify a wavenumber of $k=3$ and evaluate relative errors through the norm of (44) for each B-spline degree for meshes h0 to h4. Relative errors for this study are plotted in Figure 16 where theoretical convergence rates are demonstrated.

6.2 NASA almond

A common benchmark problem used to verify electromagnetic scattering numerical methods is the NASA almond problem as detailed in [36]. The geometry of the surface is defined through parametric expressions which are detailed in Appendix C. In the present study these expressions were used as inputs to the Math Rhino plugin developed by Rhino3DE [37] generating a NURBS representation of the almond geometry with four bicubic NURBS patches as shown in Figure 17. In addition, the software library Open CASCADE [38] was used to extract the necessary geometry data structures required to construct compatible B-spline discretizations defined over the almond surface. Div-conforming B-splines of orders $(3,2){\mkern-2.0mu\times\mkern-2.0mu}(2,3)$, $(2,1){\mkern-2.0mu\times\mkern-2.0mu}(1,2)$ and $(1,0){\mkern-2.0mu\times\mkern-2.0mu}(0,1)$ were generated with uniform h-refinement (knot insertion) applied to the initial discretization shown in Figure 17 to generate successively refined discretizations. Again, we use the notation h0, h1, h2 to indicate a mesh with no-refinement (base mesh), 1 level of h-refinement etc. and the abbreviations HH and VV to denote horizontally polarised and vertically polarised incident waves respectively. Table 2 provides further details of each B-spline discretization. For the computation of the integrals we increase the number of quadrature points in the vicinity of the two degenerate points, to increase the accuracy.

To verify our implementation we compute the RCS given by (40) at frequencies of 1.19GHz, 3GHz and 7GHz for both horizontally and vertically polarised incident waves. We use numerical RCS reference values from [39] for the 1.19GHz case, [39, 40] for the 3GHz case and [41] for the 7GHz case. In addition, we utilise experimental results for the 1.19GHz case as shown in [36]. Both [39] and [41] are based on a boundary element (method of moments) approach with the work of [40] adopting a coupled finite element/boundary element formulation.

Figure 18 illustrates RCS plots for the 1.19GHz case for each B-spline order with mesh h0. Good agreement with the numerical reference solution is visible for each order. In addition, Figure 19 demonstrates good agreement with experimental data for this frequency. In a similar manner, numerical RCS values for the 3GHz case are shown in Figure 20 where the superior accuracy of high-order discretizations is evident. Plots of the imaginary component of surface current for each order with mesh h0 are shown in Figures 21(a) to 21(c) which illustrate the smoothness in the solution obtained at higher orders.

Finally, we consider the 7GHz case where RCS plots for mesh h2 are illustrated in Figures 22 and 23 for HH and VV polarization respectively demonstrating good agreement with the numerical reference solution. At this frequency large errors were encountered for meshes h0 and h1 necessitating the use of mesh h2. Plots of the imaginary component of surface current for each order with mesh h2 are shown in Figures 24.

6.3 Integrated CAD and electromagnetic scattering analysis

We now demonstrate the ability of our approach to perform electromagnetic scattering analysis directly on CAD generated models. Figure 25 illustrates a concept model generated in Autodesk® Fusion 360™  which includes T-spline functionality capable of producing smooth, watertight surfaces. The model is composed of six bicubic NURBS surfaces consisting of 1,178 control points and 384 elements. By exporting this model as a STEP file which preserves all NURBS data structures and making use of the OpenCascade library, a compatible B-spline discretization is generated directly from this NURBS geometry model. We envisage a scenario where our implementation could be included directly with a CAD software library thereby eliminating this STEP file export procedure. The size of the bounding box for this model is given by $(\Delta x,\Delta y,\Delta z)=(82.3,93.1,27.5)$.

RCS values are computed over the $x\operatorname{-}y$ plane in which the wave is polarised in the $z$-direction. We first apply a normalised wavenumber $ka=9.31$ and apply two levels of h-refinement (denoted by $h1$ and $h2$ respectively) using compatible B-splines of order $(3,2){\mkern-2.0mu\times\mkern-2.0mu}(2,3)$ with normal $C^{0}$ continuity across patches. The discretizations $h1$ and $h2$ consist of 5,808 and 17,328 degrees of freedom respectively. Plots of the imaginary component of surface current for $h2$ are shown in Figures 26(a) and 26(b) and RCS values are plotted in Figure 27. We also compute RCS values for a higher normalised wavenumber of $ka=46.55$ in which three levels of h-refinement are applied generating a discretisation with 58,800 degrees of freedom. Surface current plots for this wavenumber are shown in Figures 28(a) and 28(b) and RCS values are plotted in Figure 29.

We use this example to demonstrate how our approach exhibits a tight link between computational design and analysis by using a common data model that provides the necessary geometry and analysis discretizations. The requirement for surface meshing is bypassed and the use of high order B-spline discretizations provides superior accuracy over conventional discretization approaches.

7 Conclusion

We have outlined an isogeometric boundary element method (method of moments) that utilises a common model to discretise both the geometry and analysis fields for electromagnetic scattering analysis. Our approach uses Non-Uniform Rational B-Splines (NURBS) to represent the surface geometry and compatible B-splines as basis for electromagnetic analysis. We have detailed the construction of compatible B-splines from a given NURBS discretization that provide a div-conforming or curl-conforming surface vector basis and described how such spline-based discretizations can be used as a basis for the electric/magnetic field integral equations. We verified our approach through the Mie series solution that provides a closed-form solution for electromagnetic scattering over a perfectly electrically conducting sphere and utilised experimental and numerical reference data for the well-known NASA almond geometry to verify radar cross section calculations. Finally, we demonstrated how our approach can be used to perform electromagnetic scattering analysis directly on geometry models generated using modern CAD software showcasing the ability of our approach to fully integrate CAD and analysis technologies.

Acknowledgements

J.A. Evans was partially supported by the Air Force Office of Scientific Research under Grant No. FA9550-14-1-0113.

Appendix A MFIE: compatible B-spline discretization

In a similar manner to the electric field integral equation, the magnetic field integral equation is first derived by substituting the expression for the total magnetic field given by

[TABLE]

into the PEC condition of

[TABLE]

to arrive at

[TABLE]

with the scattered magnetic field given by the quantity

[TABLE]

allowing (47) to be rewritten as

[TABLE]

Defining the linear operator

[TABLE]

and a forcing function $\mathbf{g}=\mathbf{n}\times\mathbf{H}^{i}$, we write the Galerkin formulation of the magnetic field integral equation as:

given $\mathbf{g}$, find $\mathbf{J}\in H^{-\frac{1}{2}}(\mathrm{curl}_{\Gamma},\Gamma)$ such that

[TABLE]

Defining finite dimensional subspaces $\mathbf{w}_{h},\mathbf{J}_{h}\in H^{-\frac{1}{2}}(\mathrm{curl}_{\Gamma},\Gamma)$ as

[TABLE]

where $\{\mathbf{N}^{\mathrm{curl}}_{A}\}_{A=1}^{N_{b}}$ is a set of curl-conforming surface vector B-spline basis functions , the system of equations for the magnetic field integral equation can be written as

[TABLE]

where, by employing the identity $\nabla\times(\phi\mathbf{v})=\nabla\phi\times\mathbf{v}+\phi\nabla\times\mathbf{v}$, applying a limiting process to the integral and noting that $\mathbf{N}_{A}^{\mathrm{div}}=-\mathbf{n}\times\mathbf{N}_{A}^{\mathrm{curl}}$,

[TABLE]

where

[TABLE]

with $\mathbf{r}:=\mathbf{y}-\mathbf{x}$ and the factor of $1/2$ arises from the limiting process. Similarly, the forcing vector components are given by

[TABLE]

As before, the vector $\mathbf{J}_{B}$ represents a vector of unknown surface current density coefficients.

Appendix B PEC sphere - additional results

Appendix C NASA almond geometry parameterization

Denoting the length of the almond geometry as $L=0.2524m$, the surface of the NASA almond geometry is defined in terms of parametric coordinates $(s,t)$ as

[TABLE]

and

[TABLE]