#!/usr/bin/env python
# coding: utf-8

# <a href='http://bempp.com' style='line-height:50px'><img src='bempp-logo.jpg' style='float:left;height:50px;margin-right:10px'> This tutorial forms part of the documentation for Bempp. Find out more at bempp.com</a>

# # Electromagnetic scattering from a screen

# ### Background

# In this tutorial, we consider the scattering of an electromagnetic wave from a perfectly conducting screen $\Gamma:=[-2,2]\times[-1,1]\times0$. The time-harmonic Maxwell equation for the electric field $\mathbf{E}$ reduces to
# 
# $$
# \nabla\times\nabla\times \mathbf{E} -k^2 \mathbf{E} = 0
# $$
# 
# in $\mathbb{R}^3\backslash\Gamma$, where $k:=2\pi/\lambda$ is the wavenumber and $\lambda$ is the wavelength. The electric field $\mathbf{E}$ is the sum of the incident field $\mathbf{E}^\text{inc}$ and the scattered field $\mathbf{E}^\text{s}$. Here, we use the incident field given by
# 
# $$
# \mathbf{E}^\text{inc}(\mathbf{x}):=\begin{bmatrix} \mathrm{e}^{\mathrm{i}kz} & 0 & 0 \end{bmatrix},
# $$
# 
# which is a wave travelling in the $z$ direction and polarised in the $x$ direction. On the screen, the tangential component $\mathbf{E}_\text{t}:=\nu\times \mathbf{E}$ must be zero. Towards infinity we impose the Silver&ndash;M&uuml;ller radiation condition
# 
# $$
# \lim_{|\mathbf{x}|\rightarrow\infty} |\mathbf{x}|\left(\nabla\times \mathbf{E}^\text{s}\times\frac{\mathbf{x}}{|\mathbf{x}|}-\mathrm{i}k\mathbf{E}^\text{s}\right) = 0,
# $$
# where $\hat{\mathbf{x}}=\mathbf{x}/|\mathbf{x}|$.
# 
# The scattered wave $\mathbf{E}^\text{s}$ has the representation
# 
# $$
# \mathbf{E}^\text{s} = -\mathcal{E}\Lambda,
# $$
# 
# where $\Lambda$ is the jump of the Neumann trace of the scattered field $\mathbf{E}^\text{s}$ across the screen. The <a href='https://bempp.com/2017/07/11/available_operators/'>Maxwell electric field potential operator</a> $\mathcal{E}$ is defined as
# 
# $$
# \mathcal{E}(\mathbf{v}):=\mathrm{i}k\int_{\Gamma}g(\mathbf{x},\mathbf{y})\mathbf{v}(\mathbf{y})\mathrm{d}\mathbf{y}-
# \frac1{\mathrm{i}k}\nabla_{\mathbf{x}}\int_{\Gamma}g(\mathbf{x},\mathbf{y})(\nabla_{\Gamma}\cdot\mathbf{v})(\mathbf{y})\mathrm{d}\mathbf{y}
# $$
# with $g(\mathbf{x},\mathbf{y}):=\frac{\mathrm{e}^{\mathrm{i}k|\mathbf{x}-\mathbf{y}|}}{4\pi|\mathbf{x}-\mathbf{y}|}$.
# 
# The associated boundary operator is denoted by $\mathsf{E}$. It is defined as average tangential trace of the electric field potential operator from both sides of the screen. The boundary integral equation is now
# 
# $$
# \mathsf{E}\Lambda = \nu\times \mathbf{E}^\text{inc}.
# $$
# The $-$ sign is missing in comparison to the representation formula since we want to satisfy the boundary conditions for the negative incident wave so that the tangential trace of the total field is zero on the screen.
# 
# More details about the mathematical background can be found in <a href='http://www.sam.math.ethz.ch/~hiptmair/Courses/CEM/BUH03.pdf' target='new'>Buffa & Hiptmair (2003)</a>.

# ### Implementation

# We start with the usual imports.

# In[1]:


import bempp.api
import numpy as np


# To avoid preconditioning issues, we assemble the operators in dense mode so that we can solve via LU decomposition later on. We will also assemble the potential operator later in dense mode.
# 
# For details of how to precondition this problem so that iterative solvers are feasible, see the <a href='https://bempp.com/2017/07/13/electric-field-integral-equation-efie/'>electric field integral equation (EFIE) tutorial</a>.

# In[2]:


bempp.api.global_parameters.assembly.boundary_operator_assembly_type = 'dense'
bempp.api.global_parameters.assembly.potential_operator_assembly_type = 'dense'


# We next define the wavenumber of the problem.

# In[3]:


wavelength = .5
k = 2 * np.pi / wavelength


# We define a structured grid in the $x$-$y$ plane.

# In[17]:


grid = bempp.api.structured_grid((-2, -1), (2, 1), (60, 30))


# We define the spaces of order 0 Raviart&ndash;Thomas div-conforming functions and order 0 N&eacute;d&eacute;lec curl-conforming functions.

# In[6]:


div_space = bempp.api.function_space(grid, "RT", 0)
curl_space = bempp.api.function_space(grid, "NC", 0)


# Next, we define the Maxwell electric field boundary operator and the identity operator. For Maxwell problems, the ``domain`` and ``range`` spaces should be div-conforming, while the ``dual_to_range`` space should be curl conforming.

# In[18]:


elec = bempp.api.operators.boundary.maxwell.electric_field(
    div_space, div_space, curl_space, k)
identity = bempp.api.operators.boundary.sparse.identity(
    div_space, div_space, curl_space)


# We create a <a href='https://bempp.com/2017/07/11/grid-functions/'>grid function</a> to represent the incident wave.

# In[22]:


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)

trace_fun = bempp.api.GridFunction(div_space, fun=tangential_trace,
                                   dual_space=curl_space)


# We use a direct LU solver to solve the system. For larger problems, <a href='https://bempp.com/2017/07/13/electric-field-integral-equation-efie/'>the problem should be preconditioned and an iterative solver used</a>.

# In[21]:


from bempp.api.linalg import lu
lambda_data = lu(elec, trace_fun)


# Now that the solution $\mathbf{\Lambda}$ is computed, we want to plot the total field. First, we define a grid of points in the $x$-$z$ plane.

# In[11]:


nx = 151
nz = 151
extent = 5
x, y, z = np.mgrid[-extent:extent:nx * 1j,
                   0:0:1j,
                   -extent:extent:nz * 1j]
points = np.vstack((x.ravel(), y.ravel(), z.ravel()))


# We now initialise the electric field potential operator.

# In[12]:


slp_pot = bempp.api.operators.potential.maxwell.electric_field(
    div_space, points, k)


# The following commands now compute the total field by first computing the scattered field from the representation formula then adding the incident field.

# In[13]:


scattered_field_data = -slp_pot * lambda_data
incident_field_data = incident_field(points)
field_data = scattered_field_data + incident_field_data


# In electromagnetic scattering it is often useful to visualise the squared electric field density. This value is computed below.

# In[14]:


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


# Finally, we can plot everything using a simple Matplotlib plot.

# In[15]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib
from matplotlib import pylab as plt
# Adjust the figure size in IPython
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0) 

plt.imshow(squared_field_density.reshape((nx, nz)).T,
           cmap='coolwarm', origin='lower',
           extent=[-extent, extent, -extent, extent])
plt.colorbar()
plt.title("Squared Electric Field Density")