Towards a Spectral Method of Moments using Computer Aided Design

We present first numerical examples of how the framework of isogeometric boundary element methods, in the context of electromagnetism also known as method of moments, can be used to achieve higher accuracies by elevation of the degree of basis functions. Our numerical examples demonstrate the computation of the electric field in the exterior domain.

1 Introduction

Spectral methods, cf. Trefethen (2000), are classes of methods for solving differential equations, closely related to classical element based methods. They share the same idea: The approximation of the solution through a series of basis functions. While classical methods choose to refine a mesh, and with this to further localise the support of each individual basis function, spectral methods employ global basis functions to approximate the solution of the problem. Since such a global basis is not always readily available, quite often any finite and boundary element method which relies solely on p-refinement, i.e., the increase of the degree of the local without mesh (h-) refinement, is referred to as a spectral element method.

In engineering applications, spectral element methods are rarely considered; the reason simply being that to fully enjoy their convergence properties, meshes with curved elements of increasing orders must be generated. This poses challenges to mesh generation and pre-processing. In contrast, classical h-refinement based mesh generators are well understood. However, with the introduction of isogeometric analysis by Hughes et al. (2005) the problem of efficient mesh generation can be avoided by the use of exact geometry mappings, allowing computations directly on CAD-generated objects.

In this document, we present numerical experiments which showcase how the isogeometric framework can be used to obtain an implementation of a spectral boundary element method. We demonstrate the implementation by the solution of electromagnetic scattering problems through p-refinement. For this, we employ a solution strategy via the electric field integral equation (Buffa and Hiptmair, 2003), which is also referred to as method of moments (MOM), and an isogeometric boundary element framework (Dölz et al., 2018b). We essentially employ p-refinement to a patchwise polynomial basis, in our case based on Bernstein polynomials. Other variants of spectral MOM have already been studied, eg. by Benoit et al. (1992) or Di Ruscio et al. (2014).

The organisation of the paper is straight forward. We first introduce basic notions of the electric field integral equation, and our p-refinement based discretization scheme, built on top of the framework of isogeometric analysis. Afterwards, we comment on the matrix assembly, followed by a discussion of our numerical examples.

Figure 1 Visualisations of the geometries. The maximum diameter of the boat is smaller than 4.5, the sphere is of radius 1.

Figure 2 Visualisations of the real part of the unknown surface current.

2 The Electric Field Integral Equation

We consider the scattering of an electromagnetic wave under the assumption of constant material coefficients μ and ϵ in Ω^c, i.e., the surroundings of a scatterer Ω, PEC boundary condition on $\mathrm{\Gamma }:=\partial \mathrm{\Omega }$ and the Silver-Müller radiation condition. Prescribing an incident wave g we arrive at

$$\begin{array}{ll}{\displaystyle }\mathbf{curl}\mathbf{curl}\mathit{e}-{\mathit{\kappa }}^{\mathrm{2}}\mathit{e}& {\displaystyle }=\mathrm{0},\mathit{\kappa }>\mathrm{0}\text{ non-resonant, in }{\mathrm{\Omega }}^{\mathrm{c}},\\ \text{(1)}& {\displaystyle }& {\displaystyle }\mathit{e}\times \mathit{n}=\mathit{g}\times \mathit{n},\end{array}$$

with wave number $\mathit{\kappa }:=\mathit{\omega }\sqrt{\mathit{ϵ}\mathit{\mu }}.$ Herein, n denotes the outward unit normal of Γ. Under these conditions, there exists a surface current j such that the scattered field can be represented by the electric field integral equation (EFIE), given by

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{e}\left(\mathit{x}\right)=-\left(\stackrel{\mathrm{̃}}{\mathfrak{V}}\mathit{j}\right)\left(\mathit{x}\right),\text{ for all }\mathit{x}\in {\mathrm{\Omega }}^{\mathrm{c}},\end{array}$$

with

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{̃}}{\mathfrak{V}}\left(\mathit{j}\right)\left(\mathit{x}\right)& {\displaystyle }:=\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}(\mathit{x},\mathit{y})\mathit{j}\left(\mathit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_y+{\displaystyle \frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}}}{\mathbf{grad}}_x\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}(\mathit{x},\mathit{y})\\ \text{(3)}& {\displaystyle }& {\displaystyle }\cdot {\mathrm{div}}_{\mathrm{\Gamma }}\left(\mathit{j}\right(\mathit{y}\left)\right)\mathrm{d}{\mathrm{\Gamma }}_y,\end{array}$$

for all $\mathit{x}\notin \mathrm{\Gamma }.$ Herein, $G_{\mathit{\kappa }}(\cdot ,\cdot )$ denotes the Green's function given by $G_{\mathit{\kappa }}(\mathit{x},\mathit{y}):=\frac{e^{i\mathit{\kappa }|\mathit{x}-\mathit{y}|}}{\mathrm{4}\mathit{\pi }|\mathit{x}-\mathit{y}|},$ cf. Buffa and Hiptmair (2003) for details.

In a continuous setting, this can be recast as the variational problem of finding an unknown surface current j in the trace of the space H(curl,Ω), such that

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{\Gamma }}{\int }\left[\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}(\mathit{x},\mathit{y})\mathit{j}\left(\mathit{x}\right)\cdot \mathit{\xi }\left(\mathit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_y+{\displaystyle \frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}}}\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}(\mathit{x},\mathit{y})\right.\\ {\displaystyle }& {\displaystyle }\left.({\mathrm{div}}_{\mathrm{\Gamma }}\circ \mathit{j})\left(\mathit{x}\right)({\mathrm{div}}_{\mathrm{\Gamma }}\circ \mathit{\xi })\left(\mathit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_y\right]\mathrm{d}{\mathrm{\Gamma }}_x=\\ \text{(4)}& {\displaystyle }& {\displaystyle }-\underset{\mathrm{\Gamma }}{\int }({\mathit{e}}_i\times \mathit{n})\cdot \mathit{\xi }\mathrm{d}\mathrm{\Gamma }\end{array}$$

holds for all test functions ξ in the same space. The corresponding trace space is often denoted by $H_{\times }^{-\mathrm{1}/\mathrm{2}}({\mathrm{div}}_{\mathrm{\Gamma }},\mathrm{\Gamma })$ and requires a divergence-conforming discretisation.

Figure 3 Numerical examples on the unit sphere. Wave number κ=1. The error refers to the maximum error obtained by comparing the values of the analytical solution with the field values through evaluation of Eq. (2) on a selection of 10 000 points on a sphere of radius 2 around the origin. Computed on a desktop PC with 16 GB RAM and an Intel i7 8700k processor.

Figure 4 Numerical examples for the ship geometry as presented in Fig. 1a. Highest p generates a discretisation with 5600 complex-valued degrees of freedom. Wave number κ=1. The error refers to the maximum error obtained by comparing the values of the analytical solution with the field values through evaluation of Eq. (2) on a selection of 10 000 points on a sphere of radius 6 around the origin. Computed on a desktop PC with 16 GB RAM and an Intel i7 8700k processor.

3 Discretisation

To solve the electric wave equation via a boundary element approach, the unknown is reduced to a vector field on Γ, often discretised by divergence-conforming elements. Implementations of such and related numerical schemes are, among others, given by Hiptmair and Kielhorn (2012), Tzoulis and Eibert (2005), or Weggler (2011). We utilise an approach based on the framework of isogeometric analysis, dealing with the special case where the B-spline basis reduces to the set of Bernstein polynomials, for order p≥0 given by

$${\displaystyle }{b}_{\mathit{\alpha },p}=\left(\begin{array}{c}n\\ \mathit{\alpha }\end{array}\right)x^{\mathit{\alpha }}(\mathrm{1}-x{)}^{p-\mathit{\alpha }},\mathrm{for}\mathrm{0}\le \mathit{\alpha }\le p.$$

To apply p-refinement, one needs sufficiently smooth (patchwise) geometry mappings, since otherwise they limit the order of ansatz functions which can be applied effectively. Thus, we choose geometry mappings given by mappings from the unit square to parts of the geometry Γ_j parametrized as tensor product mappings of rational Bernstein polynomials. Such representations can be easily extracted from NURBS (non uniform rational B-Spline) parametrisation through Bézier extraction (Borden et al., 2011).

Following Buffa et al. (2011), a suitable discretisation of the trace spaces on so-called multipatch domains $\mathrm{\Gamma }={\bigcup }_{\mathrm{0}\le j<N}{\mathrm{\Gamma }}_j$, i.e., elements in the context of spectral element methods, is provided in Buffa et al. (2018). By tensor product construction one constructs discrete spaces that are conforming w.r.t. the differential operators. As an example, if S^p denotes the Bernstein polynomials of order p≥1 on (0,1), one finds that

$${\displaystyle }{S}^p\otimes S^p\stackrel{\mathbf{curl}}{⟶}\left(\begin{array}{c}{S}^p\otimes S^{p-\mathrm{1}}\\ S^{p-\mathrm{1}}\otimes S^p\end{array}\right)\stackrel{\mathrm{div}}{⟶}{S}^{p-\mathrm{1}}\otimes S^{p-\mathrm{1}}$$

holds. This way, using the NURBS geometry mappings induced by the CAD geometries to seamlessly map between (0,1)² and patches Γ_j⊆Γ, one can define a conforming discretisation of the entire de Rham complex

$$\begin{array}{}\text{(5)}& {\displaystyle }{H}^{\mathrm{1}}\left(\mathrm{\Omega }\right)\stackrel{\mathbf{grad}}{⟶}H(\mathbf{curl},\mathrm{\Omega })\stackrel{\mathbf{curl}}{⟶}H(\mathrm{div},\mathrm{\Omega })\stackrel{\mathrm{div}}{⟶}{L}^{\mathrm{2}}\left(\mathrm{\Omega }\right),\end{array}$$

as well as its trace spaces

$$\begin{array}{}\text{(6)}& {\displaystyle }{H}^{\mathrm{1}/\mathrm{2}}\left(\mathrm{\Gamma }\right)\stackrel{{\mathbf{curl}}_{\mathrm{\Gamma }}}{⟶}{H}_{\times }^{-\mathrm{1}/\mathrm{2}}({\mathrm{div}}_{\mathrm{\Gamma }},\mathrm{\Gamma })\stackrel{{\mathrm{div}}_{\mathrm{\Gamma }}}{⟶}{H}^{-\mathrm{1}/\mathrm{2}}\left(\mathrm{\Gamma }\right),\end{array}$$

which is required for the analysis of boundary element methods. In the case of the divergence-conforming space, one must require additional normal continuity of the vector field across patch interfaces. This can be achieved through identification of the corresponding degrees of freedom (DOFs) with one another. We denote the discretisation on the boundary defined Buffa et al. (2018) by

$$\begin{array}{}\text{(7)}& {\displaystyle }{S}_p^{\mathrm{0}}\left(\mathrm{\Gamma }\right)\stackrel{{\mathbf{curl}}_{\mathrm{\Gamma }}}{⟶}{S}_p^{\mathrm{1}}\left(\mathrm{\Gamma }\right)\stackrel{{\mathrm{div}}_{\mathrm{\Gamma }}}{⟶}{S}_p^{\mathrm{2}}\left(\mathrm{\Gamma }\right).\end{array}$$

For these spaces, it is possible to show existence, uniqueness, and quasi-optimality of the solution (Dölz et al., 2018b). The discrete problem to Eq. (4) is that of finding a discrete surface current ${\mathit{j}}_h\in S_p^{\mathrm{1}}\left(\mathrm{\Gamma }\right)$ such that

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{\Gamma }}{\int }\left[\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}(\mathit{x},\mathit{y}){\mathit{j}}_h\left(\mathit{x}\right)\cdot {\mathit{\xi }}_h\left(\mathit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_{\mathit{y}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}}}\underset{\mathrm{\Gamma }}{\int }{G}_{\mathit{\kappa }}(\mathit{x},\mathit{y})({\mathbf{div}}_{\mathrm{\Gamma }}\circ {\mathit{j}}_h)\right.\\ \text{(8)}& {\displaystyle }& {\displaystyle }\left.\left(\mathit{x}\right)({\mathbf{div}}_{\mathrm{\Gamma }}\circ {\mathit{\xi }}_h)\left(\mathit{y}\right)\mathrm{d}{\mathrm{\Gamma }}_{\mathit{y}}\right]\mathrm{d}{\mathrm{\Gamma }}_{\mathit{x}}=-\underset{\mathrm{\Gamma }}{\int }({\mathit{e}}_i\times \mathit{n})\cdot {\mathit{\xi }}_h\mathrm{d}\mathrm{\Gamma }\end{array}$$

holds for all test functions ${\mathit{\xi }}_h\in S_p^{\mathrm{1}}\left(\mathrm{\Gamma }\right).$ The coefficients to represent ${\mathit{j}}_h=\sum c_j{\mathit{\xi }}_{h,j}$, where the ξ_h,j denote the basis functions of $S_p^{\mathrm{1}}\left(\mathrm{\Gamma }\right),$ can be obtained in form of the vector c given by the solution to the linear system Ac=r, where A and r can be assembled in direct analogy to Eq. (8).

4 Matrix Assembly

We apply a modified version of the superspace approach as used by e.g. Dölz et al. (2018a) to construct the dense system matrix A resulting from Eq. (4) via the representation

$${\displaystyle }\mathbf{A}={\mathbf{P}}^{\mathit{⊺}}{\mathbf{A}}^{\mathit{†}}\mathbf{P}={\mathbf{P}}^{\mathit{⊺}}\left(\sum _{i,j=\mathrm{0}}^{\mathit{\#}\text{elements}}{\mathbf{A}}_{j,i}^{\mathit{†}}\right)\mathbf{P},$$

where A^† is the dense system matrix w.r.t an elementwise polynomial, globally discontinuous basis, and P is the superspace matrix assembling the divergence-conforming basis. Herein, each matrix ${\mathbf{A}}_{j,i}^{\mathit{†}}$ is sparse, including only the interaction of the local basis on element i with that of element j. This way, the small dense system A is assembled without the need to store A^† as a whole. In the special case of only p-refinement, the matrix P incorporates only linear combinations necessary to reduce the degree in one tensor product direction for the construction of $S_p^{\mathrm{1}}\left(\mathrm{\Gamma }\right)$, and to achieve normal continuity across patch interfaces. For the solution of the arising linear system, we utilise a partially pivoted LU decomposition of the Eigen linear algebra library (Guennebaud et al., 2010).

5 Numerical Examples

To showcase the possibility of obtaining an increased accuracy through (mainly) p-refinement, we present a simple numerical example in analogy to the ones for h-refinement presented by Dölz et al. (2018a).

We excite our model setup by a Hertz-Dipole as defined by Jackson (1998, p. 411) as

$$\begin{array}{ll}{\displaystyle }{\mathit{E}}_{\mathrm{DP}}\left(\mathit{x}\right)& {\displaystyle }:=e^{i\mathit{\kappa }r}\left({\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{r}}\right(\mathit{n}\times {\mathit{p}}_{\mathrm{0}})\times \mathit{n}+({\displaystyle \frac{\mathrm{1}}{r^{\mathrm{3}}}}-{\displaystyle \frac{i\mathit{\kappa }}{r^{\mathrm{2}}}})\\ \text{(9)}& {\displaystyle }& {\displaystyle }\left(\mathrm{3}\mathit{n}\right(\mathit{n}\cdot {\mathit{p}}_{\mathrm{0}})-{\mathit{p}}_{\mathrm{0}})),\end{array}$$

with $r=\Vert \mathit{x}-{\mathit{x}}_{\mathrm{0}}\Vert$ and $\mathit{n}=(\mathit{x}-{\mathit{x}}_{\mathrm{0}})/r$. The dipole's singularity is placed inside Ω. Away from its singularity, specifically within the exterior domain, the dipole fulfils Eq. (1). Thus, by existence and uniqueness of the solution of the exterior problem Eq. (4), cf. Buffa and Hiptmair (2003), we know that Eq. (8) will approximate the surface current required to represent the field ${{\mathit{E}}_{\mathrm{DP}}}_{|{\mathrm{\Omega }}^{\mathrm{c}}},$ cf. Dölz et al. (2018b). This construction using a prescribed solution is known as a method of manufactured solutions, cf. Oberkampf and Roy (2010, Chap. 6.3).

Summarised, we solve for j in Eq. (8) with an excitation given by Eq. (9). We then evaluate Eq. (2) numerically with the approximated surface current j_h whose induced e_h(x) yields the numerical solution. By existence and uniqueness e_h is an approximation to E_DP in the exterior. We compute the error ${max}_{\mathit{x}\in V}\parallel {\mathit{E}}_{\mathrm{DP}}\left(\mathit{x}\right)-{\mathit{e}}_h\left(\mathit{x}\right){\parallel }_{C^{\mathrm{3}}}$ at points x in a set V containing points placed on a sphere enclosing the geometry, see the captions of Figs. 3 or 4 for details.

As a first example, we choose the example of a unit sphere given by 6 patches and define a dipole with ${\mathit{x}}_{\mathrm{0}}={\mathit{p}}_{\mathrm{0}}=(\mathrm{0},\mathrm{0.1},\mathrm{0.1})$ as an excitation. Although a boundary element framework has been utilised to solve the problem, no compression was applied to the systems due to their small system size, cf. Fig. 3c.

Thus the condition of the system matters little compared to the case in which compression must be applied.

On the sphere, a stable exponential rate of convergence w.r.t. p is observed. The application of p-refinement reduces the time required for matrix assembly as well as the overall system size, cf. Fig. 3, significantly.

We present another example. As geometry, we choose the boat depicted in Fig. 1 consisting of 28 Bézier patches. Placing the dipole inside with ${\mathit{x}}_{\mathrm{0}}=(\mathrm{1},\mathrm{0},\mathrm{0})$ and ${\mathit{p}}_{\mathrm{0}}=(\mathrm{0},\mathrm{0},\mathrm{0.1})$, i.e., under the “bridge”, we evaluate the electric field around the geometry. For this non-smooth and non-convex geometry, the rate of convergence, as seen in Fig. 4, is not as pronounced as before. Still, one can see an exponential convergence behaviour paired with excellent times to solution.

6 Conclusions

The adaptation of isogeometric to spectral element methods is straight forward. Through the use of Bézier extraction, one extracts piecewise smooth parametrisations of geometries. Using known constructions from isogeometric analysis, one can define a global, divergence-conforming basis that consists of tensor product Bernstein polynomials, with supports consisting of one patch each, or in the case of functions identified to achieve normal continuity, multiple patches. Application of p-refinement to this basis yields a spectral method of moments, which achieves high accuracies without the need for system compression.

These results are exceptionally promising since the reduced system size makes the application of direct solvers possible, thus circumventing the need for preconditioners, which are still a challenging topic for boundary element methods for electromagnetic problems.