#!/usr/bin/env python
# coding: utf-8
# This tutorial forms part of the documentation for Bempp. Find out more at bempp.com
# # Electric field integral equation
#
# 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 $$\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.
#
# ### Standard EFIE
# To formulate the (indirect) EFIE, we write the scattered field in the following form.
#
# $$\mathbf{E}^\text{s}=-\mathcal{E}\Lambda,$$
#
# where $\Lambda$ is an unknown tangential vector function. To find $\Lambda$, we use following boundary integral equation.
# $$\mathsf{E}\Lambda=\gamma_\mathbf{t}^\text{+}\mathbf{E}^\text{inc}.$$
#
# Here, $\gamma_\mathbf{t}^\text{+}$ is the tangential trace of a function, defined for $\mathbf{x}\in\Gamma$ by
# $$\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 $$\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
# $$\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.
# In[2]:
import bempp.api
import numpy as np
k = 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`).
# In[3]:
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.
# In[4]:
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.
# In[5]:
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.
# In[6]:
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.
# In[7]:
sol, info, iterations = bempp.api.linalg.gmres(
electric, grid_fun, return_iteration_count=True)
print("Number of iterations:", iterations)
# 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.
# In[8]:
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 we will need from it, and build a grid function representing the incident wave.
# In[9]:
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 $\mathsf{E}^2$ and $\mathsf{E}\gamma_\mathbf{t}^\text{+}\mathbf{E}^\text{inc}$.
# In[10]:
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.
# In[11]:
sol, info, iterations = bempp.api.linalg.gmres(
op, rhs, return_iteration_count=True)
print("Number of iterations:", iterations)
# 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.
# In[12]:
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.
# In[13]:
# The next command ensures that plots are shown within the IPython notebook
get_ipython().run_line_magic('matplotlib', 'inline')
# Adjust the figure size in IPython
import matplotlib
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0)
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()
# In[ ]: