# Tutorial

\[\newcommand{\vt}[1]{\boldsymbol{#1}} \newcommand{\uv}[1]{\hat{\boldsymbol{#1}}} \newcommand{\arr}[1]{\mathsf{#1}} \newcommand{\mat}[1]{\boldsymbol{\mathsf{#1}}}\]

In this tutorial we will go through the steps required for the formulation and the solution of the scattering of a time harmonic electromagnetic wave by a rectangular plate by means of the solution of the electric field integral equation.

## Building the geometry

The sibling package `CompScienceMeshes`

provides data structures and algorithms for working with simplical meshes in computational science. We will use it to create the geometry:

```
using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)
h = 0.2
Γ = meshrectangle(1.0, 1.0, h)
@show numvertices(Γ)
@show numcells(Γ)
```

```
numvertices(Γ) = 36
numcells(Γ) = 50
```

Next, we create the finite element space of Raviart-Thomas aka Rao-Wilton-Glisson functions subordinate to the triangulation `Γ`

constructed above:

`X = raviartthomas(Γ)`

The scattering problem is defined by specifying the single layer operator and the functional acting as excitation. Here, the plate is illuminated by a plane wave. The actual excitation is the tangential trace of this electric field. This trace be constructed easily by using the symbolic normal vector field `n`

defined as part of the `BEAST`

package.

```
κ = 1.0
E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
e = (n × E) × n
```

The single layer potential is also predefined by the `BEAST`

package:

`t = Maxwell3D.singlelayer(wavenumber=κ)`

It corresponds to the bilinear form

\[t(\vt{k},\vt{j}) = \frac{1}{ik} \int_{\Gamma} \int_{\Gamma'} \nabla \cdot \vt{k}(x) \nabla \cdot \vt{j}(y) \frac{e^{-ik|x-y|}}{4\pi|x-y|} dy dx - ik \int_{\Gamma} \int_{\Gamma'} \vt{k}(x) \cdot \vt{j}(y) \frac{e^{-ik|x-y|}}{4\pi|x-y|} dy dx\]

Using the `LinearForms`

package, which implements a simple form compiler for Julia (`@varform`

), the EFIE can be defined and discretised by

```
@hilbertspace j
@hilbertspace k
efie = @discretise t[k,j]==e[k] j∈X k∈X
```

Solving and computing the values of the induced current in the centers of the triangles of the mesh is now straightforward:

```
u = solve(efie)
fcr, geo = facecurrents(u,X)
```

```
dots out of 10: ..........
Converting system to Matrix...done.
LU solution of the linear system...done.
```

The resulting current distribution can be visualised by e.g. Matlab, Paraview, Plotly,...