Software frameworks for integral equations in electromagnetic scattering based on Calder\'on identities

TL;DR

This paper presents a software framework based on Calderón identities within the BEM++ library that simplifies the formulation and implementation of stable boundary integral equations for electromagnetic scattering problems.

Contribution

It introduces a robust, efficient implementation of Calderón identities and a product algebra that automates the assembly of Maxwell boundary integral equations in software.

Findings

Framework enables simple formulation of Maxwell boundary integral equations.

Implementation improves robustness and efficiency of boundary element methods.

Demonstrates practical use in deriving standard electromagnetic integral equations.

Abstract

In recent years there have been tremendous advances in the theoretical understanding of boundary integral equations for Maxwell problems. In particular, stable dual pairing of discretisation spaces have been developed that allow robust formulations of the preconditioned electric field, magnetic field and combined field integral equations. Within the BEM++ boundary element library we have developed implementations of these frameworks that allow an intuitive formulation of the typical Maxwell boundary integral formulations within a few lines of code. The basis of these developments is an efficient and robust implementations of Calder\'on identities together with a product algebra that hides and automates most technicalities involved in assembling Galerkin boundary integral equations. In this paper we demonstrate this framework and use it to derive very simple and robust software formulations of the standard preconditioned electric field, magnetic field and regularised combined field integral equations for Maxwell.

Figures

Tables

Table 1. Table 1 : Spaces to use when discretising the multitrace operator, 𝖠 𝖠 \mathsf{A} .

Matrix	Operator	Domain	Range	Dual to Range
$H_1$	Magnetic	RWG	RWG	RBC
$E_1$	Electric	BC	RWG	RBC
$E_2$	Electric	RWG	BC	SNC
$H_2$	Magnetic	BC	BC	SNC

Equations

$$\text{curl}\text{curl}{\mathbf{e}}^{\text{tot}}-k^2{\mathbf{e}}^{\text{tot}}$$

$${\mathbf{e}}^{\text{tot}}\times \mathbf{\nu }$$

$$\underset{|\mathbf{x}|\to \infty }{\lim }|\mathbf{x}|\left(\text{curl}{\mathbf{e}}^{\text{scat}}\times \frac{\mathbf{x}}{|\mathbf{x}|}-\mathrm{i}k{\mathbf{e}}^{\text{scat}}\right)$$

$${\mathbf{\gamma }}_{\text{t}}^{\pm }\mathbf{p}(\mathbf{x})$$

$${\left\{{\mathbf{\gamma }}_{\ast }\right\}}_{\Gamma }\mathbf{f}$$

$${\mathbf{L}}_{\text{t}}^2(\Gamma ):=\{\mathbf{u}\in {\mathbf{L}}^2(\Gamma ):\mathbf{u}\cdot \mathbf{\nu }=0\},$$

$${\mathbf{H}}_{\times }^{1/2}(\Gamma ):={\mathbf{\gamma }}_{\text{t}}^{\text{–}}({\mathbf{H}}^1({\Omega }^{\text{–}}))=\left\{{\mathbf{\gamma }}_{\text{t}}^{\text{–}}\mathbf{u}:\mathbf{u}\in {\mathbf{H}}^1({\Omega }^{\text{–}})\right\}.$$

$${⟨\mathbf{a},\mathbf{b}⟩}_{\tau }$$

$${\text{div}}_{\Gamma }({\mathbf{\gamma }}_{\text{t}}^{\pm }\mathbf{u}):={\mathbf{\gamma }}_{\mathbf{\nu }}^{\pm }(\text{curl}\mathbf{u}),$$

$${\mathbf{\gamma }}_{\mathbf{\nu }}^{\pm }\mathbf{p}:=\underset{{\Omega }^{\pm }\ni {\mathbf{x}}^{\prime }\to \mathbf{x}\in \Gamma }{\lim }\mathbf{p}({\mathbf{x}}^{\prime })\cdot \mathbf{\nu }(\mathbf{x}).$$

$${\text{div}}_{\Gamma }\mathbf{u}:=\{\begin{array}{ll}{\text{div}}_j{\mathbf{u}}^j& \text{on }{\Gamma }^j\\ ({\mathbf{u}}^j\cdot {\mathbf{\nu }}^{ij}+{\mathbf{u}}^i\cdot {\mathbf{\nu }}^{ji}){\delta }_{ij}& \text{on }\overline{{\Gamma }^j}\cap \overline{{\Gamma }^i}\end{array},$$

$${\mathbf{H}}_{\times }^{-1/2}({\text{div}}_{\Gamma },\Gamma ):=\{\mathbf{\mu }\in {\mathbf{H}}_{\times }^{-1/2}(\Gamma ):{\text{div}}_{\Gamma }\mathbf{\mu }\in H^{-1/2}(\Gamma )\}.$$

$${\text{curl}}_{\Gamma }(\mathbf{u}):={\text{div}}_{\Gamma }(\mathbf{u}\times \mathbf{\nu }),$$

$${\mathbf{H}}_{\times }^{-1/2}({\text{curl}}_{\Gamma },\Gamma ):=\{\mathbf{\nu }\times \mathbf{\mu }:\mathbf{\mu }\in {\mathbf{H}}_{\times }^{-1/2}({\text{div}}_{\Gamma },\Gamma )\}.$$

$${\mathbf{\gamma }}_{\text{t}}^{\pm }$$

$$\text{and }{\mathbf{\gamma }}_{\text{N}}^{\pm }$$

$${\mathbf{RT}}_i(\mathbf{r}):=\{\begin{array}{ll}\frac{1}{2A_i^{\text{+}}}(\mathbf{r}-{\mathbf{p}}_i^{\text{+}})& \text{if }\mathbf{r}\in T_i^{\text{+}}\\ -\frac{1}{2A_i^{\text{–}}}(\mathbf{r}-{\mathbf{p}}_i^{\text{–}})& \text{if }\mathbf{r}\in T_i^{\text{–}}\\ 0& \text{otherwise}\end{array},$$

$${\mathbf{NC}}_i(\mathbf{r}):=\mathbf{\nu }\times {\mathbf{RT}}_i(\mathbf{r}).$$

$${\mathbf{RWG}}_i(\mathbf{r}):=l_i{\mathbf{RT}}_i(\mathbf{r})=\{\begin{array}{ll}\frac{l_i}{2A_i^{\text{+}}}(\mathbf{r}-{\mathbf{p}}_i^{\text{+}})& \text{if }\mathbf{r}\in T_i^{\text{+}}\\ -\frac{l_i}{2A_i^{\text{–}}}(\mathbf{r}-{\mathbf{p}}_i^{\text{–}})& \text{if }\mathbf{r}\in T_i^{\text{–}}\\ 0& \text{otherwise}\end{array},$$

$${\mathbf{SNC}}_i(\mathbf{r}):=l_i{\mathbf{NC}}_i(\mathbf{r})=\mathbf{\nu }\times {\mathbf{RWG}}_i(\mathbf{r}).$$

$$\underset{\mathbf{u}\in \text{RWG}}{\inf }\underset{\mathbf{v}\in \text{BC}}{\sup }\frac{⟨\mathbf{u},\mathbf{v}{⟩}_{\tau }}{\Vert \mathbf{u}{\Vert }_{-\frac{1}{2},\text{div}}\Vert \mathbf{v}{\Vert }_{-\frac{1}{2},\text{div}}}>\frac{1}{C}.$$

$${\mathbf{RBC}}_i(\mathbf{r}):=\mathbf{\nu }\times {\mathbf{BC}}_i(\mathbf{r}).$$

$$\mathcal{E}(\mathbf{p})(\mathbf{x})$$

$$\mathcal{H}(\mathbf{p})(\mathbf{x})$$

$$\mathbf{e}(\mathbf{x})=-\mathcal{H}([{\mathbf{\gamma }}_{\text{t}}{]}_{\Gamma }\mathbf{e})(\mathbf{x})-\mathcal{E}([{\mathbf{\gamma }}_{\text{N}}{]}_{\Gamma }\mathbf{e})(\mathbf{x}).$$

$$\mathsf{E}$$

$${\left\{{\mathbf{\gamma }}_{\text{N}}\right\}}_{\Gamma }\mathcal{E}$$

$${\left[{\mathbf{\gamma }}_{\text{t}}\right]}_{\Gamma }\mathcal{E}$$

$${\mathbf{\gamma }}_{\text{t}}^{\text{+}}\mathcal{E}$$

$${\mathbf{\gamma }}_{\text{t}}^{\text{+}}\mathcal{H}$$

$${\mathbf{\gamma }}_{\text{t}}^{\text{–}}\mathcal{E}$$

$${\mathbf{\gamma }}_{\text{t}}^{\text{–}}\mathcal{H}$$

$$\mathsf{A}:=\left[\begin{array}{cc}\mathsf{H}& \mathsf{E}\\ -\mathsf{E}& \mathsf{H}\end{array}\right].$$

Full text

††journal: CAMWA

1 Introduction

The numerical simulation of electromagnetic wave scattering poses significant theoretical and computational challenges. Much effort in recent years has gone into the development of fast and robust boundary integral equation formulations to simulate a range of phenomena from the design and performance of antennas to radar scattering from large metallic objects.

While there have been a range of important theoretical advances in recent years for the development of robust preconditioned boundary integral formulations for Maxwell, the computational implementation remains a challenge. At University College London, as part of the BEM++ project [1] we have developed a number of easy to use Python-based open-source tools to explore and solve Maxwell problems based on preconditioned electric field (EFIE), magnetic field (MFIE) and combined field (CFIE) integral equation formulations. In this paper, we give an overview of these developments, present simple example codes and discuss the underlying implementation. The goal is to make advanced integral equation solvers for Maxwell available to non-specialised users without requiring significant knowledge in the design and implementation of these methods.

In particular, we consider the following formulation of electromagnetic scattering from a perfectly conducting object. Let $\Omega^{\textup{--}}\subset\mathbb{R}^{3}$ be a bounded domain with boundary $\Gamma$ and denote by $\Omega^{\textup{+}}=\mathbb{R}^{3}\backslash\overline{\Omega^{\textup{--}}}$ its complement. Let ${\bm{\nu}}$ be the exterior normal vector on $\Gamma$ pointing into $\Omega^{\textup{+}}$ as shown in Figure 1.

Denote by ${\bm{e}}^{\textup{inc}}$ an incident field. We are looking for the solution ${\bm{e}}^{\textup{tot}}={\bm{e}}^{\textup{inc}}+{\bm{e}}^{\textup{scat}}$ of the exterior scattering problem, satisfying

[TABLE]

Here, $k=\omega\sqrt{\epsilon_{0}\mu_{0}}$ denotes the wavenumber of the problem, with $\omega$ denoting the frequency and $\epsilon_{0}$ and $\mu_{0}$ the electric permeability and magnetic permittivity in vacuum. Frequently, the incident field is a plane wave given by ${\bm{e}}^{\textup{inc}}={\bm{p}}e^{\mathrm{i}k{\bm{x}}\cdot{\bm{d}}}$, where ${\bm{p}}$ is a non-zero vector representing the polarisation of the wave and ${\bm{d}}$ is a unit vector perpendicular to ${\bm{p}}$ that gives the direction of the plane wave.

In Section 2, we give a short overview of the main concepts of the BEM++ boundary element library with an emphasis on its operator and grid function algebra.

When discretising these equations, the finite dimensional function spaces must be chosen carefully in order to give well-conditioned systems of equations. In Section 3, we present an overview of the function spaces for Maxwell boundary integral operators, and the various discretisations of these spaces. The emphasis is on the presentation of stable dual pairings between so-called Rao–Wilton–Glisson (RWG) [2] spaces, and Buffa–Christiansen (BC) spaces [3].

In Section 4, we present the standard electric and magnetic field boundary integral operators and the resulting Calderón projector. The concept of the Calderón projector is fundamental to this paper. We describe the function spaces used in assembling the Calderón projector and numerically demonstrate its properties using small BEM++ code snippets.

In Sections 5, 6 and 7, we then demonstrate the Calderón preconditioned EFIE, MFIE and regularised CFIE and their BEM++ implementations. Interesting non-trivial domains will be used to demonstrate and compare their properties.

The use of stable RWG/BC pairings, on which this paper is based, is well known in the electromagnetics community (see e.g. [4, 5, 6]). The emphasis in this paper is on a simple high-level software representation that is mathematically correct and robust, but hides the underlying complexities to provide a simple framework in which to solve electromagnetic scattering problems.

2 A Maxwell introduction to BEM++

BEM++ (www.bempp.org) is an open-source library for the Galerkin discretisation and solution of boundary integral equations in electrostatics, acoustics and computational electromagnetics. The library supports dense discretisation and hierarchical matrix assembly. The interface of the library is written in Python with fast kernel assembly routines implemented in C++.

In this section we give a brief overview of BEM++ using the example of Maxwell boundary integral operators. We skip mathematical details: these will be discussed in the subsequent sections. A particular focus in this section is on the operator concept of BEM++, which will allow us to elegantly describe the various types of integral equation in what follows.

BEM++ has a range of features to define or import triangular surface grids consisting of flat elements. To define a function space over a grid we use the function_space command. For example, a simple space consisting of RWG functions (see Section 3 for details) is defined by

⬇

import bempp.api

grid = …

rwg_space = bempp.api.function_space(grid, 'RWG', 0)

The third parameter in the function_space command is the degree of the space. For Maxwell only lowest order functions are supported.

With the definition of a function space we can now discretise functions over a given space into a GridFunction object. The following code computes a discretisation of the tangential trace of a plane wave Maxwell solution ${\bm{u}}({\bm{x}})={\bm{p}}e^{\mathrm{i}k{\bm{d}}\cdot{\bm{x}}}$, where ${\bm{p}}$ is the polarisation vector.

⬇

import numpy as np

p = np.array([…])

d = np.array([…])

k = …

def plane_wave(x):

**return** np.exp(1j * np.dot(d, x)) * p

def tangential_trace(x, n, domain_index, result):

result[:] = np.cross(plane_wave(x), n, axis=0)

grid_fun = bempp.api.GridFunction(

    space=rwg_space, fun=tangential_trace)

The tangential_trace function defines the tangential trace of the plane wave in dependence of a given normal vector n, which is passed in during the discretisation phase. The GridFunction object by default performs an ${\bm{L}}^{2}$ projection onto the space rwg_space. Oblique projections with a different test space are also supported by specifying the dual_space argument in the constructor of GridFunction.

We now define an operator acting on grid functions. An electric field integral operator is defined on the space ${\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)$ of surface-div-conforming functions. For the discretisation we need test functions from the space of surface-curl-conforming functions ${\bm{H}}_{\times}^{-1/2}(\textup{curl}_{\Gamma},\Gamma)$. Here, we choose RBC functions, which are explained in Section 3.

⬇

from bempp.api.operators.boundary import \

    maxwell

rbc_space = bempp.api.function_space(grid, 'RBC', 0)

op = maxwell.electric_field(

    rwg_space, rwg_space, rbc_space, k)

Each boundary operator in BEM++ takes three parameters, a domain space, a range space, and a dual space to the range space. For a Galerkin discretisation the range space is not required. However, we use it to implement a Galerkin based operator algebra. This algebra allows us for example to write

⬇

result = op * grid_fun

The result is a grid function, which is defined over the range space of the operator. Internally, this is implemented by multiplying the weak form of the operator with the coefficient vector of the grid function, and then solving with the mass matrix between the dual space and the range space to map the result back into a coefficient vector of the range space. This mass matrix solve in practice is only performed if a method requires the coefficient vector of the result grid function, so that the above code snippet is even well-defined if there is no stable mapping from the dual space to the range space. However, then only a representation of the grid function in the dual space is available.

By extending this mechanism to operators, the concept of products of operators is also defined in BEM++. Given two operators op1 and op2 in BEM++ with compatible spaces, their product is a new operator that is implicitly defined through first the application of op2, followed by a mass matrix solve to map the result into the domain space of op1, followed by the application of op1.

The high-level concept of operators and grid functions extends everywhere into the library. For example, to solve the indirect electric field integral equation (EFIE) with the above plane wave incident field we can write

⬇

sol = bempp.api.linalg.lu(op, grid_fun)

to solve the associated EFIE problem using a dense LU decomposition. The object sol is a grid function that lives in the domain space of op and is a grid function that satisfies grid_fun == op * sol up to numerical errors. Naturally, for larger problems we are interested in iterative solvers and corresponding interfaces to GMRES and CG are also provided.

We will use this operator concept extensively in this paper to formulate the different types of Maxwell integral equations. More details of the abstract implementation of operator algebras in BEM++ can be found in [7].

3 Function spaces

In this section we review the necessary function space definitions for boundary integral formulations of the electromagnetic scattering problem (1). In what follows, we assume that $\Gamma$ consists of a finite number of smooth faces that meet at non-degenerate edges. This is a reasonable assumption from the point of view of boundary element methods, as we will consider discrete problems on meshed surfaces.

The fundamental difficulty lies in the correct description of trace spaces on $\Gamma$ for solutions of (1) [8, 9, 10, 11, 12]. In this section, we summarise without proof some of the results of these papers, as they form the foundations for the operator definitions and stable implementations presented in later sections. Throughout, we will present small code snippets to demonstrate how to instantiate the corresponding spaces in BEM++. More details and references can also be found in the overview paper [13] and in [14].

Throughout this paper, we adopt the convention of using bold and non-bold letters for spaces of vector- and scalar-valued functions, respectively. Whenever we write the subscript loc it is understood to only be used for the unbounded domain $\Omega^{+}$. For definitions of Sobolev spaces of scalar-valued functions, we refer to [15, chapter 3].

Let op be one of curl, $\textup{{curl}}^{2}$ or div. Denote by ${\bm{H}}_{\text{loc}}(\textbf{op},\Omega^{\pm})$ the function space defined as ${\bm{H}}_{\text{loc}}(\textbf{op},\Omega^{\pm}):=\{{\bm{u}}\in{\bm{L}}^{2}_{\text{loc}}(\Omega^{\pm}):\textbf{op}\,{\bm{u}}\in{\bm{L}}^{2}_{\text{loc}}(\Omega^{\pm})\}$ (for div, the final ${\bm{L}}^{2}_{\text{loc}}(\Omega^{\pm})$ should be replaced by the scalar space $L^{2}_{\text{loc}}(\Omega^{\pm})$). To define the necessary function spaces on the surface $\Gamma$, we first define the tangential (${}_{\text{t}}$) and Neumann (${}_{\text{N}}$) traces on $\Gamma$. These are defined, for ${\bm{p}}\in{\bm{H}}_{\text{loc}}(\textup{{curl}},\Omega^{\pm})$ and ${\bm{q}}\in{\bm{H}}_{\text{loc}}(\textup{{curl}}^{2},\Omega^{\pm})$, by

[TABLE]

where the superscripts ${}^{\textup{--}}$ and ${}^{\textup{+}}$ denote the interior and exterior traces, respectively. Note that in our definition ${\bm{\gamma}}_{\textup{N}}^{\pm}$ contains an additional factor of $\mathrm{i}$, which does not appear in [13]. The interpretation is that if we normalise the magnetic permittivity and electric permeability to $1$, this definition of ${\bm{\gamma}}_{\textup{N}}^{\pm}$ is the tangential trace of the magnetic field data.

In what follows we need the average $\left\{\cdot\right\}_{\Gamma}$, and jump, $\left[\cdot\right]_{\Gamma}$ of these traces, defined as

[TABLE]

Let ${\bm{L}}^{2}_{\textup{t}}(\Gamma)$, defined by

[TABLE]

be the space of square integrable tangential functions. We define the tangential trace space, ${\bm{H}}_{\times}^{1/2}(\Gamma)$, as in [13, definition 1] by

[TABLE]

The dual of this space with respect to the antisymmetric product,

[TABLE]

is denoted by ${\bm{H}}_{\times}^{-1/2}(\Gamma)$.

Due to the assumption that $\Gamma$ consists of a finite number of smooth faces, we may let $\Gamma=\bigcup_{j=0}^{m}\Gamma^{j}$, where $\Gamma_{0}$, …, $\Gamma_{m}$ are smooth. We define the scalar surface divergence of ${\bm{\gamma}}_{\textup{t}}^{\pm}{\bm{u}}$, for ${\bm{u}}\in{\bm{H}}_{\text{loc}}(\textup{{curl}},\Omega^{\pm})$, by

[TABLE]

where the normal trace, ${\bm{\gamma}}_{{\bm{\nu}}}^{\pm}$, is defined for ${\bm{p}}\in{\bm{H}}_{\text{loc}}(\textup{{div}},\Omega^{\pm})$ by

[TABLE]

For a function ${\bm{u}}\in{\bm{\gamma}}_{\textup{t}}^{\textup{--}}({\bm{C}}^{\infty}(\overline{\Omega^{\textup{--}}}))$, this can be calculated using

[TABLE]

where ${\bm{u}}^{j}$ is the restriction of ${\bm{u}}$ to the face $\Gamma^{j}$, ${\bm{\nu}}^{ij}$ is the outward pointing tangential normal to $\Gamma^{i}$ restricted to the edge $\overline{\Gamma^{i}}\cap\overline{\Gamma^{j}}$, $\textup{div}_{j}$ is the two dimensional divergence computed on the face $\Gamma^{i}$, and $\delta_{ij}$ is the Dirac delta distribution with support on the edge $\overline{\Gamma^{i}}\cap\overline{\Gamma^{j}}$. By density, this definition can be extended to ${\bm{u}}\in{\bm{H}}_{\times}^{-1/2}(\Gamma)$. We now define the space of surface-div-conforming functions by

[TABLE]

The scalar surface curl may be defined, for ${\bm{u}}\in{\bm{H}}^{-1/2}_{\times}(\Gamma)$, by

[TABLE]

and the space of surface-curl-conforming functions by

[TABLE]

By restricting their domains as in the following theorem, we obtain continuous and surjective traces.

Theorem 1**.**

The traces

[TABLE]

are continuous and surjective.

Proof.

See [11, Theorem 4.1] and [13, Lemma 3]. ∎

The antisymmetric dual form defined in (6) is intimately connected with the space ${\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)$. In [11, Lemma 5.6] it is shown that the space ${\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)$ is self-dual with respect to this antisymmetric dual form. Another interpretation of this dual form is as the standard ${\bm{L}}^{2}$ dual between the spaces ${\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)$ and ${\bm{H}}_{\times}^{-1/2}(\textup{curl}_{\Gamma},\Gamma)$ since ${\bm{\psi}}\in{\bm{H}}_{\times}^{-1/2}(\textup{curl}_{\Gamma},\Gamma)$ if and only if ${\bm{\psi}}={\bm{\nu}}\times{\bm{\xi}}$ for some ${\bm{\xi}}\in{\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)$.

Commonly used discretisations of these two spaces are Raviart-Thomas (RT) div-conforming [16] and Nédélec (NC) curl-conforming [17] basis functions.

For the $i$th edge in a mesh, between two triangles $T_{i}^{\text{+}}$ and $T_{i}^{\text{--}}$, the order 0 RT basis function is defined by

[TABLE]

where $A_{i}^{\text{+}}$ and $A_{i}^{\text{--}}$ are the areas of $T_{i}^{\text{+}}$ and $T_{i}^{\text{--}}$, and ${\bm{p}}_{i}^{\text{--}}$ and ${\bm{p}}_{i}^{\text{+}}$ are the corners of $T_{i}^{\text{+}}$ and $T_{i}^{\text{--}}$ not on the shared edge, as shown in Figure 2. For the same edge, the order 0 NC basis function may be defined by

[TABLE]

Example RT and NC order 0 basis functions are shown in Figure 3. In BEM++, RT and NC spaces may be created with the following lines of Python.

⬇

rt_space = bempp.api.function_space(grid, "RT", 0)

nc_space = bempp.api.function_space(grid, "NC", 0)

The RT basis functions are closely related to the Rao–Wilton–Glisson (RWG) basis functions presented in [2]. These are defined by

[TABLE]

where $l_{i}$ is the length of the shared edge, and all other terms are as above. We define the scaled curl-conforming dual basis functions of the RWG functions as

[TABLE]

In BEM++, RWC and SNC spaces may be created with the following lines of Python.

⬇

rwg_space = bempp.api.function_space(grid, "RWG", 0)

snc_space = bempp.api.function_space(grid, "SNC", 0)

An overview of other bases that can be used to discretise ${\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)$ and ${\bm{H}}_{\times}^{-1/2}(\textup{curl}_{\Gamma},\Gamma)$, as well as other spaces, can be found in [18].

3.1 Buffa–Christiansen spaces

The RT basis functions have a subspace that is quasi-orthogonal to the space of curl-conforming Nédélec functions [19, section 3.1]. Due to this, the antisymmetric bilinear form, as defined in (6), on the discrete RT space is not inf-sup stable. The motivation for Buffa–Christiansen (BC) basis functions is to find a space of functions that are div-conforming but behave like curl-conforming functions, as this will recover inf-sup stability.

To define such a basis, we first take the space of div-conforming RWG functions on a barycentric refinement of the original grid. Given the $i$th edge of the original coarse grid, the associated BC function $\mathbf{BC}_{i}$ is now defined as a linear combination of these RWG functions on the barycentric refinement such that this linear combination is approximately tangentially (and therefore curl-) conforming across the $i$th edge. A typical BC basis function is shown in Figure 5. The full definition of BC functions is given in [3] in which also the following inf-sup stability result is shown. This result implies the stability of the Gram matrix between BC functions and RWG functions with respect to the dual pairing in (6).

[TABLE]

For a BC basis function, $\mathbf{BC}_{i}$, we may also define the rotated Buffa–Christiansen (RBC) basis function, in an analogous way to (16), by

[TABLE]

In BEM++, BC and RBC spaces may be created with the following lines of Python.

⬇

bc_space = bempp.api.function_space(grid, "BC", 0)

rbc_space = bempp.api.function_space(grid, "RBC", 0)

We can use BEM++ to compare the stability of dual pairings of RWG spaces with SNC and RBC spaces. The code in Figure 6 computes the condition number of mass matrices on a regular sphere grid generated from the RWG/SNC pairing (id1) and the RWG/RBC pairing (id2).

For the condition number of id1 the code computes a value of $7.7\times 10^{17}$ and for the condition number of id2 it computes a value of $3.60$. We note that in the definitions of the spaces in Figure 6 we have used the identifiers B-RWG and B-SNC instead of RWG and SNC. The reason is that the RBC spaces are defined over barycentric refinements of the grid. So we need to tell also the other space definitions to internally use barycentric refinements of the grid (even though the actual spaces live on the coarse grid), which is done by prepending B- in the definitions.

4 Operators

BEM++ provides the magnetic and electric field domain potential and boundary operators.

First, we define the electric and magnetic potential operators (see [13]), $\mathcal{E},\mathcal{H}:{\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)\to{\bm{H}}_{\text{loc}}(\textup{{curl}}^{2},\Omega^{\textup{+}}\cup\Omega^{\textup{--}})$, by

[TABLE]

where $\displaystyle G({\bm{x}},{\bm{y}})=\frac{e^{\mathrm{i}k|{\bm{x}}-{\bm{y}}|}}{4\pi|{\bm{x}}-{\bm{y}}|}$ is the Green’s function of a three-dimensional Helmholtz problem.

The definition of the electric potential operator, $\mathcal{E}$, used here differs from that used in [13] by a factor of $\mathrm{i}$, corresponding to the modified definition of ${\bm{\gamma}}_{\textup{N}}^{\pm}$.

With the electric and magnetic field operators we obtain the following representation formula.

Theorem 2**.**

If ${\bm{e}}\in{\bm{H}}_{\textup{loc}}(\textup{{curl}}^{2},\Omega^{\textup{--}}\cup\Omega^{\textup{+}})$ is a solution of (1a), then

[TABLE]

Proof.

See [13, section 4]. ∎

Once the jumps of the traces of the solution are known or estimated on $\Gamma$, the representation formula (23) can be used to find the solution at points in $\Omega^{\pm}$. In BEM++ the electric and magnetic domain potential operators are available to evaluate this representation formula.

⬇

from bempp.api.operators.potential import maxwell

eval_points = …

e_pot = maxwell.electric_field(

    magnetic_space, eval_points, k)

m_pot = maxwell.magnetic_field(

    electric_space, eval_points, k)

field = - e_pot * magnetic_jump - m_pot * electric_jump

The domain is represented by a $3\times N$ array of $N$ points away from the boundary $\Gamma$. By default the potential operators are assembled using $\mathcal{H}$-matrices and then for given trace data evaluated by a fast $\mathcal{H}$-matrix/vector product [20].

Taking traces of the electric and magnetic field potential operators we arrive at the electric boundary operator, $\mathsf{E}:{\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)\to{\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)$, and the magnetic boundary operator, $\mathsf{H}:{\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)\to{\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)$. These are defined by

[TABLE]

Additionally, we define the identity operator, $\mathsf{Id}$, that maps every function to itself. Because of the symmetry between electric and magnetic fields, the average Neumann traces can be written in terms of $\mathsf{E}$ and $\mathsf{H}$ as follows:

[TABLE]

The following jump conditions can be derived [13, Theorem 7].

[TABLE]

Combining (24),(25) and (26) gives

[TABLE]

for the exterior traces and

[TABLE]

for the interior traces.

For given domain, range_ and dual_to_range spaces the electric and magnetic field boundary operators in BEM++ are defined by

⬇

from bempp.api.operators.boundary import maxwell

efie = maxwell.electric_field(

    domain, range_, dual_to_range, k)

mfie = maxwell.magnetic_field(

    domain, range_, dual_to_range, k)

To obtain standard RWG discretisations of these operators we choose RWG as domain and range space, and SNC as dual to range space. It is important to note that internally we use the formulations of the weak forms of these operators given in [13], which are based on the antisymmetric dual form $\langle\cdot,\cdot\rangle_{\tau}$ and not on the standard ${\bm{L}}^{2}$ dual form. Hence, the dual spaces are the non-rotated spaces RWG and BC for the discretisation, and when SNC is passed as dual to range space, it is internally interpreted as an RWG space for the generation of the discrete weak form. The rotated spaces SNC and RBC play a role in the operator algebra used to form the mass matrices for the $\langle\cdot,\cdot\rangle_{\tau}$ dual form by means of the standard ${\bm{L}}^{2}$ dual form.

4.1 The Calderón projector and its discretisation

Using the boundary operators in the previous section, we define the multitrace operator $\mathsf{A}$ by

[TABLE]

We then define the exterior Calderón projector, $\mathcal{C}^{\textup{+}}$, as follows.

[TABLE]

It is well known [13] that, given any arbitrary trace data $({\bm{a}},{\bm{b}})\in{\bm{H}}_{\times}^{-1/2}(\textup{div}_{\Gamma},\Gamma)^{2}$ the product

[TABLE]

defines a compatible pair of Cauchy data, and

[TABLE]

from which we obtain $\left[\mathcal{C}^{\textup{+}}\right]^{2}=\mathcal{C}^{\textup{+}}$. Using this identity and the representation $\mathcal{C}^{\textup{+}}=\tfrac{1}{2}\mathsf{Id}-\mathsf{A}$ we obtain

[TABLE]

This relationship is crucial for preconditioning numerical methods based on the Calderón projector and any discretisation scheme should preserve this property.

Denote by

[TABLE]

the discretisation of the operator $\mathsf{A}$. Here, $E_{1}$ and $E_{2}$ are discretisations of electric field operators and $H_{1}$ and $H_{2}$ are discretisations of magnetic field operators.

We want to numerically satisfy the relationship in (33). Hence, we require that

[TABLE]

where $M_{1}$ and $M_{2}$ are the corresponding mass matrices between the dual spaces and range spaces of the operators in the first line, and in the second line, respectively.

In order to satisfy the above relationship it is crucial that $M_{1}$ and $M_{2}$ are well-conditioned mass matrices. A choice of spaces for the operators that achieves this goal, is shown in Table 1. These choices of spaces lead to all mass matrices in the discretisation of $\mathsf{A}^{2}$ being the invertible RWG-RBC or BC-SNC pairings. The choice of spaces in Table 1 is based on representing the tangential trace with an RWG space and the Neumann trace with a BC space. Alternatively, one could use BC for the electric component and RWG for the magnetic component. This would lead to a discretisation in which $E_{1}$ and $E_{2}$ are swapped, and $H_{1}$ and $H_{2}$ are swapped.

Using the BEM++ library, the stable multitrace operator may be created using the following lines of Python.

⬇

from bempp.api.operators.boundary import maxwell

multitrace = maxwell.multitrace_operator(grid, k)

We may then create the exterior Calderón projector with the following lines.

⬇

from bempp.api.operators.boundary import sparse

identity = sparse.multitrace_identity(

grid, spaces="maxwell")

calderon = 0.5 * identity - multitrace

A complete example for a computation with the Calderón projector is given in Figure 7. It takes as electric and magnetic trace data the tangential trace of the constant vector $(1,0,0)$. This is obviously not a pair of compatible Cauchy data for Maxwell. It then applies the Calderón projector to obtain the numerically compatible pair of Cauchy data traces_1 and then applies the Calderón projector one more time to obtain the pair of traces traces_2. In a stable implementation, both, traces_1 and traces_2 should agree up to discretisation error, and indeed, the error electric_error in the electric component is $9.8\times 10^{-3}$ and the error magnetic_error in the magnetic component is $7.4\times 10^{-3}$.

It is important to note that we need to choose different discretisation spaces for the electric and magnetic trace. For the electric trace we choose RWG basis functions, and for the magnetic trace we choose BC basis functions.

In Figure 8, we show the electric and magnetic trace data obtained from taking the tangential trace of the vector $(1,0,0)$ (upper plot), and the result of applying the Calderón projector to the incompatible pair of Cauchy data (bottom plot).

4.2 Implementational Details

The discrete multitrace operator $A$ consists of the two magnetic field operator discretisations $H_{1}$ and $H_{2}$, and the two electric field operator discretisations $E_{1}$ and $E_{2}$. In practice, we only create two operators $\tilde{E}$, and $\tilde{H}$, using RWG basis functions on the barycentrically refined grid. Let $\tilde{M}$ be the $L^{2}$ mass matrix associated with this RWG space, defined by

[TABLE]

where $\mathbf{RWG}_{i}$ is the $i$th RWG basis function on the barycentrically refined grid. Let now $M_{RWG}$ be the mass matrix obtained from trial functions in the RWG space on the barycentrically refined grid, and test functions from the original RWG space on the coarse grid. We correspondingly define the mass matrix $M_{BC}$ with test functions from the BC space on the original coarse grid. The operators $E_{i}$, and $H_{i}$, $i=1,2$, are now given as

[TABLE]

In [4], a similar construction of the matrices is suggested. The difference is that there permutation matrices that represent the basis functions on the coarse mesh in terms of basis functions on the barycentric refinement are stated explicitly.

The implicit construction here has the advantage that it is independent of the particular space. All that is needed is the ability to construct mass matrices, which is often already available. A potential performance pitfall is the application of the mass matrix inverse of $\tilde{M}$ for each matrix vector product with $E_{i}$ or $E_{j}$. We automatically precompute the LU decomposition of $\tilde{M}$. Even for fairly large meshes this is done in a few seconds.

5 Electric field integral equation (EFIE)

We now turn our attention to the EFIE. The EFIE is widely used in applications for low-frequency scattering from closed and open surfaces. In its strong form the indirect EFIE to find the scattered field ${\bm{e}}^{\textup{scat}}$ in (1) can be written as (see [13])

[TABLE]

Correspondingly, the direct EFIE can be written as

[TABLE]

The direct EFIE is derived from the first line of the exterior Calderón projector. The unknown ${\bm{\pi}}$ has a physical meaning as the Neumann trace of the scattered field. Once, the solution ${\bm{\lambda}}$ for the indirect EFIE or ${\bm{\pi}}$ for the direct EFIE is available the solution ${\bm{e}}^{\textup{scat}}$ can be computed using

[TABLE]

for the indirect EFIE, and

[TABLE]

for the direct EFIE.

By itself the EFIE is ill-conditioned, making necessary either direct solvers or efficient preconditioners. The Calderón identities described in Section 4.1 provide an efficient preconditioning strategy. From the top-left block of $\mathsf{A}^{2}=\frac{1}{4}\mathsf{Id}$ it follows that $\mathsf{E}^{2}=-\tfrac{1}{4}\mathsf{Id}+\mathsf{H}^{2}$. The eigenvalues of $\mathsf{E}$ accumulate at [math] and $\infty$ making discretisations of this operator highly ill-conditioned. However, the operator $\mathsf{H}$ is compact on smooth surfaces [21, Section 5.5]. Hence, the eigenvalues of $\mathsf{E}^{2}$ accumulate at $-\frac{1}{4}$. Implementations of this self-regularising property depend on the ability to perform the operator product $\mathsf{E}^{2}$ in a stable way. In Section 4.1, we saw how stable operator products are implemented, based on Buffa–Christiansen bases, and this property will be used here to implement the preconditioning strategy. Further details on the Calderón preconditioner for the EFIE can also be found in [4].

The numerical implementation of the EFIE is typically based on discretising the operator $\mathsf{E}$ in (35) and (36) by pairs of RWG/SNC trial and dual spaces. For the direct EFIE, this implies that a discretisation of the Calderón projector is used, in which the tangential trace data is represented with BC basis functions and the Neumann trace data is represented with RWG basis functions.

In Figure 9 we show the BEM++ implementation of the Calderón preconditioned indirect EFIE, based on the stable formulation of the multitrace operator $\mathsf{A}$. The EFIE operator is $E_{2}$ from the discretisation of $\mathsf{A}$, which maps from the RWG to the BC space. The preconditioning operator is $E_{1}$, which maps from the BC space to the RWG space. Correspondingly, the right-hand side incident wave is discretised using BC functions to be compatible with the preconditioning operator. The solution sol lives in the RWG space, which is the domain space of $E_{2}$. The GMRES routine has the additional parameter use_strong_form=True. This enables mass matrix preconditioning and further reduces the number of iterations. At the end we discretise a domain potential operator over the RWG space of the solution function sol. This evaluates the potential at arbitrary points in the domain.

We emphasise that BEM++ automatically only assembles a single electric field operator on the barycentric refinement of the grid. The magnetic components of the multitrace operator are not assembled as they are not needed for the solution of the problem.

Taking $\Omega$ as unit sphere, the incident wave $\displaystyle{\bm{e}}^{\textup{inc}}=\left[e^{\mathrm{i}kz},0,0\right]$, and $k=2$, we discretise the EFIE on a series of triangular grids with different levels of refinement. The number of GMRES iterations required to solve the linear system arising from EFIE and the Calderón preconditioned EFIE applied to this example problem are shown in Figure 10. It can be seen that while the number of iterations required to solve the EFIE rises quickly as the grid is refined, the number of iterations required to solve the preconditioned EFIE remains below 10.

Figures 11 and 12 show electromagnetic waves scattering off more interesting obstacles. Both were computed using the indirect Calderón preconditioned EFIE.

The left plot of Figure 11 shows the wave ${\bm{e}}^{\textup{inc}}=[0,0,e^{\mathrm{i}kx}]$, with $k=20\pi$, scattering off the NASA almond benchmarking shape, as defined in [22]. This was discretised on a grid with 2442 edges. The right plot shows the corresponding convergence curve of the GMRES residual for the preconditioned EFIE.

Figure 12 shows the wave ${\bm{e}}^{\textup{inc}}={\bm{p}}e^{\mathrm{i}k{\bm{d}}\cdot{\bm{x}}}$, with ${\bm{p}}=[-1,2,0]$, ${\bm{d}}=[\tfrac{2}{\sqrt{5}},\tfrac{1}{\sqrt{5}},0]$, and $k=5$, scattering off a level 1 Menger sponge, and the corresponding convergence curve. The Menger sponge was discretised using a grid with 4680 edges.

6 Magnetic field integral equation (MFIE)

The MFIE can be represented in BEM++ as easily as the EFIE. The strong form of the indirect MFIE is (see [13])

[TABLE]

It is a valid formulation on closed domains. Its advantage compared to the EFIE is that on smooth domains, it is a compact perturbation of the identity operator and therefore well suited to iterative solvers. However, the robust implementation of the MFIE also on non-smooth domains requires the use of RWG trial spaces and RBC test spaces (see [5]).

The MFIE can be easily implemented as shown in the code snippet in Figure 13.

In Figure 14, we demonstrate the rate of convergence of the MFIE for solving the NASA almond and the Menger sponge examples from Section 5. Both converge nicely without additional preconditioning. While the rate of convergence is slower than in the preconditioned EFIE case, only one boundary operator needs to be applied for each iteration step, while the preconditioned EFIE requires two applications.

7 Combined field integral equation (CFIE)

While the EFIE and MFIE are efficient for low-frequency Maxwell problems, they lead to break-down close to interior resonances. The CFIE is immune to breakdown at resonances and is therefore particularly suitable for high-frequency scattering problems. Here, we focus on the direct CFIE and the stable version of it derived in [23]. In its strong form, it is given as

[TABLE]

The CFIE is a linear combination of the direct MFIE, obtained from the second row of the exterior Calderón projector, and a regularised direct EFIE, obtained by multiplying the first row of the exterior Calderón projector with a regularisation operator $\mathsf{R}$. Frequently, the EFIE component is multiplied with a complex scalar. This is not necessary here, as in our implementation the electric field operator itself is already scaled with $\mathrm{i}$.

We use the RWG space for the unknown Neumann trace ${\bm{\pi}}$. Hence, we swap $E_{1}$ and $E_{2}$, and $H_{1}$ and $H_{2}$ in the discretisation of the Calderón projector in Section 4.1.

It follows that we discretise $\mathsf{H}$ on the left-hand side of (39) with $H_{1}$. Moreover, for $\mathsf{E}$ on the left-hand side we choose the discrete operator $E_{2}$. The operator $E_{2}$ maps from RWG into BC, while $H_{1}$ maps from RWG into RWG. We therefore require that a discretisation of $\mathsf{R}$ maps from the BC space to the RWG space. We could for example choose the operator $E_{1}$. But this operator is not injective at interior electric eigenvalues. The solution is to choose $E_{1}$ based on the wavenumber $\mathrm{i}k$, instead of $k$ (see [23]). On the right-hand side we choose for $\mathsf{H}$ the discretisation $H_{2}$ and for $\mathsf{E}$ the discretisation $E_{1}$ to stay compatible with the corresponding direct EFIE and direct MFIE formulation. We can easily implement this in the framework of BEM++ with the code snippet in Figure 15.

In Figure 16, we demonstrate our CFIE implementation for a simple scattering problem on the unit sphere for varying wavenumbers. While the CFIE shows a bounded and only slowly growing number of GMRES iterations, the number of iterations for the EFIE and MFIE grow strongly in the neighborhood of an interior resonance close to $k=12.5$. In Figure 17, we show the convergence of GMRES for the CFIE applied to the NASA almond and Menger sponge examples from Sections 5 and 6.

Another construction of the CFIE based on the use of BC spaces was presented in [24]. In particular, the treatment of the MFIE component in that paper differs from the proposed formulation in this section.

8 Concluding remarks

We have only discussed the Maxwell solver aspects of BEM++. The goal is to make available advanced techniques for the solution of electromagnetic scattering problems while hiding complex details associated with robust boundary integral equation formulations for Maxwell. We have not discussed here the solution of Maxwell transmission problems. But the proposed framework fits this very well as a we need for each medium the corresponding multitrace operator $\mathsf{A}$. Work on single-trace and multi-trace transmission formulations is ongoing.

We have not included timing benchmarks in this paper since the focus is on the demonstration of the software framework and not the performance comparison of different solvers. BEM++ internally uses its own thread and MPI parallelised $\mathcal{H}$-matrix implementation for the discretisation of integral operators, which is very efficient for low-frequency problems, but also performs well in the medium frequency range.

Finally, we want to conclude this paper with an example adapted to the theme of this special issue. Figure 18 shows the electromagnetic field around perfectly conducting birthday cake (not for eating!) computed with BEM++.

All examples in this paper were computed with version 3.1 of the BEM++ library, and will be made available in the form of IPython Notebooks on www.bempp.org.

Acknowledgements

The authors would like to thank Peter Monk for his continuing support of the BEM++ project. Indeed, he was the first external user of BEM++ when the project just started, and over the years he has constantly given helpful feedback and bug reports.

We would furthermore like to thank Dave Hewett (UCL) and Andrea Moiola (Reading) for helpful discussions regarding trace spaces of boundary integral operators for Maxwell.

The work of the second author was supported by Engineering and Physical Sciences Research Council grant EP/K03829X/1.