#!/usr/bin/env python
# coding: utf-8
# This tutorial forms part of the documentation for Bempp. Find out more at bempp.com
# # 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).
#
# ## Background
#
# 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.
#
# ### MFIE
# The (indirect) MFIE uses the representation formula
#
# $$\mathbf{E}^\text{s}=-\mathcal{H}\Lambda,$$
#
# and the following boundary integral equation.
#
# $$\left(\mathsf{H}-\tfrac12\mathsf{Id}\right)\Lambda=\gamma_\mathbf{t}^\text{+}\mathbf{E}^\text{inc}.$$
#
# 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.
# ## Implementation
# 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,
fun=tangential_trace,
dual_space=rbc_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)