Magnetic field integral equation

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


In this tutorial, we use consider an incident wave $$\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 $\mathbf{E}^\text{s}$ be the scattered field and look to solve $$ \begin{align} \textbf{curl}\,\textbf{curl}\,\mathbf{E}-k^2\mathbf{E}&=0\quad\text{in }\Omega^\text{+},\\ \mathbf{E}\times\nu&=0\quad\text{on }\Gamma,\\ \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, \end{align} $$

where $\mathbf{E}=\mathbf{E}^\text{s}+\mathbf{E}^\text{inc}$ is the total electric field. This is the same problem as solved in the EFIE example.


The (indirect) MFIE uses the representation formula


and the following boundary integral equation.


Here, $\gamma_\mathbf{t}^\text{+}$ is the tangential trace of a function, as defined in the EFIE example.

In Cools et al (2011), it was suggested that the robust implementation of the MFIE on non-smooth domains requires the use of the stable space pairings, as described in the multitrace example script.


First, we do the usual imports, set the wavenumber, and define the incident wave, as in the EFIE example.

In [1]:
import bempp.api
import numpy as np

k = 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 and operator we will need from it, and build a grid function representing the incident wave.

In [2]:
multitrace = bempp.api.operators.boundary.maxwell.multitrace_operator(
    grid, k)
identity = bempp.api.operators.boundary.sparse.multitrace_identity(
    grid, spaces="maxwell")
calderon = .5 * identity - multitrace

rwg_space = multitrace.domain_spaces[0]
rbc_space = multitrace.dual_to_range_spaces[0]

rhs = bempp.api.GridFunction(rwg_space,
op = -calderon[0,0]

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

In [3]:
sol, info, iterations = bempp.api.linalg.gmres(
    op, rhs, return_iteration_count=True)
print("Number of iterations:", iterations)
Number of iterations: 20