#!/usr/bin/env python
# coding: utf-8
# 
This tutorial forms part of the documentation for Bempp. Find out more at bempp.com
# # Combined field integral equation
#
# This tutorial shows how to solve the combined field integral equation (CFIE) for exterior scattering problems, as described in section 7 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 in the EFIE and MFIE examples.
#
# ### The CFIE
# Both the MFIE and the EFIE are efficient for low-frequency problems, but at higher frequencies they are unstable close to interior resonances of the domain. The CFIE is immune to these resonances, so is well suited to high-frequency problems.
#
# The CFIE is a linear combination of the MFIE and EFIE. Its strong form is as follows:
#
# $$-\mathsf{R}\mathsf{E}\Lambda + (\tfrac12\mathsf{Id}+\mathsf{H})\Lambda
# =-\mathsf{R}(\tfrac12\mathsf{Id}+\mathsf{H})\gamma_\mathbf{t}^\text{+}\mathbf{E}^\text{inc}-\mathsf{E}\gamma_\mathbf{t}^\text{+}\mathbf{E}^\text{inc},$$
#
# where $\mathsf{R}$ is a regularisatior operator. In this example, we use $\mathsf{R}=\mathsf{E}_{\mathrm{i}k}$, the electric field boundary operator with wavenumber $\mathrm{i}k$. More details of the CFIE can be found in Contopanagos et al (2002).
#
# As we did with the EFIE and MFIE, we need to carefully choose the spaces used to discretise the CFIE to ensure that all the products of operators are stable.
# ## Implementation
# First, we do the usual imports, set the wavenumber, and define the incident wave, as in the EFIE and MFIE examples.
# 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 operators with wavenumbers $k$ and $\mathrm{i}k$, extract the spaces and operators we will need from them, and build a grid function representing the incident wave.
# In[2]:
multitrace = \
bempp.api.operators.boundary.maxwell.multitrace_operator(
grid, k)
multitrace_scaled = \
bempp.api.operators.boundary.maxwell.multitrace_operator(
grid, 1j*k)
identity = bempp.api.operators.boundary.sparse.multitrace_identity(
grid, spaces='maxwell')
rwg_space = multitrace.domain_spaces[0]
snc_space = multitrace.dual_to_range_spaces[1]
bc_space = multitrace.domain_spaces[1]
rbc_space = multitrace.dual_to_range_spaces[0]
calderon = .5 * identity + multitrace
grid_fun = bempp.api.GridFunction(
bc_space, fun=tangential_trace,
dual_space=snc_space)
R = multitrace_scaled[0, 1]
E1 = multitrace[0, 1]
E2 = -multitrace[1, 0]
mfie1 = calderon[0, 0]
mfie2 = calderon[1, 1]
# We use the operators to form the two sides of the equation.
# In[3]:
op = -R * E2 + mfie1
rhs = -R * mfie2 * grid_fun - E1 * grid_fun
# Next, we define the incendent field and its tangential trace.
# In[4]:
co, info, its = bempp.api.linalg.gmres(
op, rhs, use_strong_form=True,
return_iteration_count=True)
print("Number of iterations:", its)
# By running this example for a range of wavenumbers, we can see that the CFIE is indeed immune to resonances. In the following plot, the CFIE is represented by black diamonds, the MFIE by blue squares, and the Calderón preconditioned EFIE by red circles.
# 