Electric field integral equation (EFIE)

Download this tutorial as a Jupyter notebook

This tutorial shows how to solve the electric field integral equation (EFIE) for exterior scattering problems, as described in section 5 of Scroggs et al (2017).

Background

In this tutorial, we use consider an incident wave

\displaystyle \mathbf{E}^\text{inc}(\mathbf{x})=\left(\begin{array}{c}\mathrm{e}^{\mathrm{i}kz}\\0\\0\end{array}\right)

scattering off the unit sphere.

We let \displaystyle \mathbf{E}^\text{s} be the scattered field and look to solve

\displaystyle   \textbf{curl}\,\textbf{curl}\,\mathbf{E}-k^2\mathbf{E}=0\quad\text{in }\Omega^\text{+},\\[2mm]  \mathbf{E}\times\nu=0\quad\text{on }\Gamma,\\[2mm]  \lim_{|\mathbf{x}|\to\infty}\left(\textbf{curl}\,\mathbf{E}^\text{s}\times\frac{\mathbf{x}}{|\mathbf{x}|}-\mathrm{i}k\mathbf{E}^\text{s}\right)=0,

where \displaystyle \mathbf{E}=\mathbf{E}^\text{s}+\mathbf{E}^\text{inc} is the total electric field.

Standard EFIE

To formulate the (indirect) EFIE, we write the scattered field in the following form.

\displaystyle \mathbf{E}^\text{s}=-\mathcal{E}\Lambda,

where \displaystyle \Lambda is an unknown tangential vector function. To find \displaystyle \Lambda, we use following boundary integral equation.

\displaystyle \mathsf{E}\Lambda=\gamma_\mathbf{t}^\text{+}\mathbf{E}^\text{inc}.

Here, \displaystyle \gamma_\mathbf{t}^\text{+} is the tangential trace of a function, defined for \displaystyle \mathbf{x}\in\Gamma by

\displaystyle \gamma_\mathbf{t}^\text{+}\mathbf{A}(\mathbf{x}) := \lim_{\Omega^\text{+}\ni\mathbf{x'}\to\mathbf{x}}\mathbf{A}(\mathbf{x}')\times\nu(\mathbf{x}).

Calderón preconditioned EFIE

The boundary integral equation for the EFIE is ill-conditioned, and so will be infeasible to solve for larger problems. Using properties of the multitrace operator, we can show that

\displaystyle \mathsf{E}^2=-\tfrac14\mathsf{Id}+\mathsf{H}.

This is a compact perturbation of the identity, and so this will lead to a well conditioned system.

The boundary integral equation for the Calderón preconditioned EFIE is therefore

\displaystyle \mathsf{E}^2\Lambda=\mathsf{E}\gamma_\mathbf{t}^\text{+}\mathbf{E}^\text{inc}.

As mentioned in the multitrace operator tutorial, the spaces used to discretise this must be chosen carefully to ensure that a stable discretisation is achieved.

Implementation

First, we do the usual imports and set the wavenumber.

import bempp.api
import numpy as np

= 3

Next, we define the grid. In the paper, we use the sphere, plus the Nasa almond (bempp.api.shapes.almond) and a level 1 Menger sponge (bempp.api.shapes.menger_sponge).

grid = bempp.api.shapes.sphere(h=0.1)

We will first solve the non-preconditioned EFIE. For this, we define the spaces of Raviart–Thomas (RT) and Nédélec (NC) functions.

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

Next, we define the incendent field and its tangential trace.

def incident_field(x):
    return np.array([np.exp(1j*k*x[2]), 0.*x[2], 0.*x[2]])

def tangential_trace(x, n, domain_index, result):
    result[:] = np.cross(incident_field(x), n, axis=0)

grid_fun = bempp.api.GridFunction(rt_space,
                                  fun=tangential_trace,
                                  dual_space=nc_space)

We define the electric field operator, using RT functions for the domain and range spaces and NC functions for the dual space.

electric = bempp.api.operators.boundary.maxwell.electric_field(
    rt_space, rt_space, nc_space, k)

Finally, we solve the discretisation of the problem and print the number of iterations.

sol, info, iterations = bempp.api.linalg.gmres(
    electric, grid_fun, return_iteration_count=True)

print("Number of iterations:", iterations)

Number of iterations: 522

As expected, the number of iterations taken to solve the non-preconditioned system is high.

Calderón preconditioned EFIE

To solve the preconditioned EFIE, we begin by importing Bempp and Numpy and defining the wavenumber and incident wave as above.

import bempp.api
import numpy as np

= 3

grid = bempp.api.shapes.sphere(h=0.1)

def incident_field(x):
    return np.array([np.exp(1j*k*x[2]), 0.*x[2], 0.*x[2]])

def tangential_trace(x, n, domain_index, result):
    result[:] = np.cross(incident_field(x), n, axis=0)

We define the multitrace operator, extract the spaces we will need from it, and build a grid function representing the incident wave.

multitrace = bempp.api.operators.boundary.maxwell.multitrace_operator(
    grid, k)
bc_space = multitrace.range_spaces[1]
snc_space = multitrace.dual_to_range_spaces[1]

grid_fun = bempp.api.GridFunction(bc_space,
                                  fun=tangential_trace,
                                  dual_space=snc_space)

We extract the electric field operators E1 and E2 from the multitrace operator, then form the products \displaystyle \mathsf{E}^2 and \displaystyle \mathsf{E}\gamma_\mathbf{t}^\text{+}\mathbf{E}^\text{inc}.

E2 = -multitrace[1,0]
E1 = multitrace[0,1]
op = E1 * E2
rhs = E1 * grid_fun

Next, we solve the discrete system and print the number of iterations.

sol, info, iterations = bempp.api.linalg.gmres(
    op, rhs, return_iteration_count=True)
print("Number of iterations:", iterations)
Number of iterations: 13

As expected, the preconditioned system requires a much lower number of iterations.

To plot a slice of the solution, we define a grid of points and use the representation formula to evaluate the squared electric field density at these points.

x_p, y_p, z_p = np.mgrid[-5:5:300j, 0:0:1j, -5:5:300j]

points = np.vstack((x_p.ravel(), y_p.ravel(), z_p.ravel()))

efie_pot = bempp.api.operators.potential.maxwell.electric_field(
    sol.space, points, k)
plot_me = incident_field(points) - efie_pot * sol

plot_me = np.real(np.sum(plot_me * plot_me.conj(), axis=0))

for i,p in enumerate(points.T):
    if np.linalg.norm(p) <= 1:
        plot_me[i] = None

Finally, we plot the slice of the solution.

from matplotlib import pyplot as plt

plt.imshow(plot_me.reshape((300,300)).T, origin='lower',
           extent=[-5,5,-5,5], vmin=0, vmax=4, cmap='coolwarm')
plt.colorbar()
plt.title("Plot at y=0")
plt.xlabel("x")
plt.ylabel("z")
plt.show()