In this tutorial, we consider the scattering of an electromagnetic wave from a perfectly conducting screen . The time-harmonic Maxwell equation for the electric field reduces to
in , where is the wavenumber and is the wavelength. The electric field is the sum of the incident field and the scattered field . Here, we use the incident field given by
which is a wave travelling in the direction and polarised in the direction. On the screen, the tangential component must be zero. Towards infinity we impose the Silver–Müller radiation condition
The scattered wave has the representation
where is the jump of the Neumann trace of the scattered field across the screen. The Maxwell electric field potential operator is defined as
The associated boundary operator is denoted by . 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
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 Buffa & Hiptmair (2003).
We start with the usual imports.
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 electric field integral equation (EFIE) tutorial.
bempp.api.global_parameters.assembly.potential_operator_assembly_type = ‘dense’
We next to define the wavenumber of the problem.
k = 2 * np.pi / wavelength
We define a structured grid in the – plane.
We define the spaces of order 0 Raviart–Thomas div-conforming functions and order 0 Nédélec curl-conforming functions.
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.
div_space, div_space, curl_space, k)
identity = bempp.api.operators.boundary.sparse.identity(
div_space, div_space, curl_space)
We create a grid function to represent the incident wave.
return np.array([np.exp(1j * k * x), 0. * x, 0. * x])
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,
We use a direct LU solver to solve the system. For larger problems, the problem should be preconditioned and an iterative solver used.
lambda_data = lu(elec, trace_fun)
Now that the solution is computed, we want to plot the total field. First, we define a grid of points in the – plane.
nz = 151
extent = 5
x, y, z = np.mgrid[–extent:extent:nx * 1j,
–extent:extent:nz * 1j]
points = np.vstack((x.ravel(), y.ravel(), z.ravel()))
We now initialise the electric field potential operator.
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.
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.
Finally, we can plot everything using a simple Matplotlib plot.
from matplotlib import pylab as plt
extent=[–extent, extent, –extent, extent])
plt.title(“Squared Electric Field Density”)