The multitrace operator and Calderón projector

Download this tutorial as a Jupyter notebook

This tutorial shows how to create and use a Calderón projectors and multitrace operators for Maxwell problems, as described in section 4 of Scroggs et al (2017).

Background

For Maxwell problems, we define the multitrace operator \displaystyle \mathsf{A} by

\displaystyle \mathsf{A}=\begin{bmatrix}\mathsf{H}&\mathsf{E}\\-\mathsf{E}&\mathsf{H}\end{bmatrix},

where \displaystyle \mathsf{H} and \displaystyle \mathsf{E} are the magnetic and electric field boundary operators.

The interior Calderón projector \displaystyle \mathcal{C}^\text{--} is defined by

\displaystyle \mathcal{C}^\text{--}=\tfrac12\mathsf{Id}+\mathsf{A},

and the exterior Calderón projector \displaystyle \mathcal{C}^\text{+} is defined by

\displaystyle \mathcal{C}^\text{+}=\tfrac12\mathsf{Id}-\mathsf{A}.

If \displaystyle \mathbf{E} and \displaystyle \mathbf{H} are valid Cauchy data for a exterior Maxwell problem, then

\displaystyle \mathcal{C}^\text{+}\begin{bmatrix}\mathbf{E}\\\mathbf{H}\end{bmatrix}=\begin{bmatrix}\mathbf{E}\\\mathbf{H}\end{bmatrix}.

It follows from the property \displaystyle \left[\mathcal{C}^\text{+}\right]^2=\mathcal{C}^\text{+} that for any tangential functions \displaystyle \mathbf{A} and \displaystyle \mathbf{B}, the functions \displaystyle \mathbf{E} and \displaystyle \mathbf{H} defined by

\displaystyle \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 \displaystyle \mathsf{B} and \displaystyle \mathsf{C} be two boundary operators. If the product \displaystyle \mathsf{BC} were discretised, the result would \displaystyle B_hM^{-1}C_h, where \displaystyle M is the identity matrix between the range and dual_to_range spaces of \displaystyle \mathsf{C}.

If standard Raviart–Thomas (RT) and Nédélec (NC) basic functions are used to discretise \displaystyle \mathsf{C}, the matrix \displaystyle M will be numerically singular. To avoid this problem, a stable pairing of RT and Buffa–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 \displaystyle \mathsf{A}^2 are stable. This is discussed in more detail in section 4 of Scroggs et al (2017).

Implementation

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

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

= 2

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

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.

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

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

calderon = 0.5 * identity - multitrace

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

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.

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 \displaystyle \left[\mathcal{C}^\text{+}\right]^2=\mathcal{C}^\text{+}. To verify this, we plot the relative difference between the functions.

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)
Electric error is 0.00983125547609
Magnetic error is 0.00740308354962

As expected, the errors are small.