#!/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>

# # The multitrace operator and Calderón projector
# 
# This tutorial shows how to create and use a Calderón projectors and multitrace operators for Maxwell problems, as described in section 4 of <a href='https://bempp.com/documentation/example-scripts-maxwell/'>Scroggs <em>et al</em> (2017)</a>.
# 
# ## Background
# For Maxwell problems, we define the multitrace operator $\mathsf{A}$ by
# 
# $$\mathsf{A}=\begin{bmatrix}\mathsf{H}&\mathsf{E}\\-\mathsf{E}&\mathsf{H}\end{bmatrix},$$
# 
# where $\mathsf{H}$ and $\mathsf{E}$ are the <a href='https://bempp.com/2017/07/11/available_operators/'>magnetic and electric field boundary operators</a>.
# 
# The interior Calder&oacute;n projector $\mathcal{C}^\text{--}$ is defined by
# 
# $$\mathcal{C}^\text{--}=\tfrac12\mathsf{Id}+\mathsf{A},$$
# 
# and the exterior Calder&oacute;n projector $\mathcal{C}^\text{+}$ is defined by
# 
# $$\mathcal{C}^\text{+}=\tfrac12\mathsf{Id}-\mathsf{A}.$$
# 
# If $\mathbf{E}$ and $\mathbf{H}$ are valid Cauchy data for a exterior Maxwell problem, then
# 
# $$\mathcal{C}^\text{+}\begin{bmatrix}\mathbf{E}\\\mathbf{H}\end{bmatrix}=\begin{bmatrix}\mathbf{E}\\\mathbf{H}\end{bmatrix}.$$
# 
# It follows from the property $\left[\mathcal{C}^\text{+}\right]^2=\mathcal{C}^\text{+}$ that for any tangential functions $\mathbf{A}$ and $\mathbf{B}$, the functions $\mathbf{E}$ and $\mathbf{H}$ defined by
# 
# $$\begin{bmatrix}\mathbf{E}\\\mathbf{H}\end{bmatrix}=\mathcal{C}^\text{+}\begin{bmatrix}\mathbf{A}\\\mathbf{B}\end{bmatrix}$$
# 
# are valid Cauchy data.
# 
# ### Discretisation
# Let $\mathsf{B}$ and $\mathsf{C}$ be two boundary operators. If the product $\mathsf{BC}$ were discretised, the result would $B_hM^{-1}C_h$, where $M$ is the identity matrix between the ``range`` and ``dual_to_range`` spaces of $\mathsf{C}$.
# 
# If standard Raviart&ndash;Thomas (RT) and N&eacute;d&eacute;lec (NC) basic functions are used to discretise $\mathsf{C}$, the matrix $M$ will be numerically singular. To avoid this problem, a stable pairing of RT and Buffa&ndash;Christiansen (BC) functions should be used.
# 
# In Bempp, the function ``bempp.api.operators.boundary.maxwell.multitrace_operator`` will automatically use the correct spaces so that all terms in the product $\mathsf{A}^2$ are stable. This is discussed in more detail in section 4 of <a href='https://bempp.com/documentation/example-scripts-maxwell/'>Scroggs <em>et al</em> (2017)</a>.
# 
# ## Implementation

# First, we import everything we need and set the wavenumber.

# In[1]:


import bempp.api
from bempp.api.operators.boundary import maxwell
from bempp.api.operators.boundary import sparse
import numpy as np

k = 2


# We make the meshed cube on which we will define our operators.

# In[2]:


grid = bempp.api.shapes.cube(h=0.1)


# Next, we make the multitrace operator and identity on the grid. Here, we must rell the multitrace identity to use vector ``"maxwell"`` spaces instead of the default scalar spaces used for Laplace and Helmholtz.

# In[3]:


multitrace = maxwell.multitrace_operator(grid, k)
identity = sparse.multitrace_identity(grid, spaces='maxwell')


# We define the exterior Calderón projector, $\mathcal{C}^\text{+} = \tfrac12\mathsf{Id}-\mathsf{A}$.

# In[4]:


calderon = 0.5 * identity - multitrace


# We define grid functions using a Python function. These functions will not be valid Cauchy data.

# In[5]:


def tangential_trace(x, n, domain_index, result):
    result[:] = np.cross(np.array([1,0,0]),n)

electric_trace = bempp.api.GridFunction(
    space=calderon.domain_spaces[0],
    fun=tangential_trace,
    dual_space=calderon.dual_to_range_spaces[0])

magnetic_trace = bempp.api.GridFunction(
    space=calderon.domain_spaces[1],
    fun=tangential_trace,
    dual_space=calderon.dual_to_range_spaces[1])


# We now apply the Calderón projector to the functions twice.

# In[ ]:


traces_1 = calderon * [electric_trace, magnetic_trace]
traces_2 = calderon * traces_1


# As explained above, the functions in ``traces_1`` should be valid Cauchy data. Additionally, ``traces_2`` should be equal to ``traces_1`` because $\left[\mathcal{C}^\text{+}\right]^2=\mathcal{C}^\text{+}$. To verify this, we plot the relative difference between the functions.

# In[ ]:


electric_error = (traces_2[0] - traces_1[0]).l2_norm() / traces_1[0].l2_norm()
magnetic_error = (traces_2[1] - traces_1[1]).l2_norm() / traces_1[1].l2_norm()
print("Electric error is ", electric_error)
print("Magnetic error is ", magnetic_error)


# As expected, the errors are small.

# In[ ]:


traces_1[0].plot()


# In[ ]: