Helmi-FEM 0.3.0

Finite element method for solving the electromagnetic Helmholtz equation.

Helmi_FEM

An open-source implementation of the Helmholtz equation for finite element analysis of electromagnetic wave propagation and scattering.

Based on scikit-fem for finite element assembly and on SciPy, NumPy, and matplotlib for solving, processing and plotting.

Installation

pip install Helmi-FEM

Usage

import helmi

See the example workflow below and additional examples in the subfolder.

Features

• Supports complex-valued parameters
• TE and TM polarization
• Dirichlet and third-order boundary conditions (Neumann, freespace, absorbing)
• Near-field to far-field transformation

Workflow

1. Create or load a mesh of the simulation domain

Use the mesh constructors of skfem.mesh for simple shapes such as rectangles, circles, or L-shapes. Use gmsh or other software for complex shapes.

https://scikit-fem.readthedocs.io/en/latest/api.html#module-skfem.mesh

http://gmsh.info/

import skfem
import numpy as np
from helmi import Helmholtz

# create a rectangular mesh with skfem:
x_pts = np.linspace(0, 100, 101)
y_pts = np.linspace(-5, 5, 21)
mesh = skfem.MeshTri.init_tensor(x_pts, y_pts)
mesh = mesh.with_subdomains({'air': lambda x: x[0] < 50,
                             'plastic': lambda x: x[0] >= 50})
mesh = mesh.with_boundaries({'bound_xmin': lambda x: np.isclose(x[0], x_pts[0]),
                             'bound_xmax': lambda x: np.isclose(x[0], x_pts[-1]),
                             'bound_ymin': lambda x: np.isclose(x[1], y_pts[0]),
                             'bound_ymax': lambda x: np.isclose(x[1], y_pts[-1])})

# alternatively, load a mesh from file:
mesh = skfem.Mesh.load('mymesh.msh')

2. Choose a finite element from skfem.element.

https://scikit-fem.readthedocs.io/en/latest/api.html#module-skfem.element

element = skfem.ElementTriP2()
fem = Helmholtz(mesh, element)

3. Assemble the domains and the boundaries

The Helmholtz equation is implemented as a general second-order partial differential equation:

$-\alpha (\frac{{\partial }^2\mathrm{\Phi }}{\partial x^2}+\frac{{\partial }^2\mathrm{\Phi }}{\partial y^2})+\beta \mathrm{\Phi }=f$.

The electromagnetic wave propagation can be formulated in two different ways, depending on the polarization of interest:

1. TM polarization: $\mathrm{\Phi }=E_Z$, $\alpha =1/{\mu }_r$, $\beta =-k_0^2{ϵ}_r$, and $f=-j{k}_0{Z}_0{J}_Z$
2. TE polarization: $\mathrm{\Phi }=H_Z$, $\alpha =1/{ϵ}_r$, $\beta =-k_0^2{\mu }_r$, and $f=\frac{1}{{ϵ}_r}(\frac{\partial J_Y}{\partial x}-\frac{\partial J_X}{\partial y})$

with the normalized free-space vacuum wave number $k_0=\frac{2\pi }{{\lambda }_0}=\frac{2\pi f}{c_0}$. Here, $c_0$ is propagation velocity (speed of light) in vacuum, expressed in mesh units per second.

For $f=10GHz$ and a mesh unit of $a=1mm$, one gets $k_0\approx \frac{2\pi 10GHz}{3e8m/s/1mm}\approx 0.209$.

For a mesh with two different subdomains labeled 'air' and 'plastic', the domains are assembled with their respective parameters:

k0 = 0.5
eps_air = 1
mu_air = 1
eps_plastic = 2 - 0.1j
mu_plastic = 1
fem.assemble_subdomains(alpha={'air': 1 / mu_air, 
                               'plastic': 1 / mu_plastic}, 
                        beta={'air': -1 * k0 ** 2 * eps_air, 
                              'plastic': -1 * k0 ** 2 * eps_plastic}, 
                        f={'air': 0, 
                           'plastic': 0})

Similarly, the boundary conditions are defined. In this example, the upper and lower boundaries, labeled as 'bound_ymax' and 'bound_ymin', are supposed to be perfectly conducting (metallic) waveguide walls. The corresponding essential (Dirichlet) boundary condition is $E_Z=0$:

fem.assemble_boundaries_dirichlet(value={'bound_ymin': 0, 
                                         'bound_ymax': 0})

The boundaries on the left and right, labeled 'bound_xmin' and 'bound_xmax', are supposed to be waveguide interfaces that are artificially extended to behave like infinite waveguides for the fundamental propagation mode. This can be formulated with the third-order boundary condition: $\alpha (\frac{\partial \mathrm{\Phi }}{\partial x}\hat{x}+\frac{\partial \mathrm{\Phi }}{\partial y}\hat{y})\hat{n}+\gamma \mathrm{\Phi }=q$. Here, $\hat{x}$, $\hat{y}$, and $\hat{n}$ are unit vectors in x and y direction, and outward normal to the boundary, respectively.

The desired waveguide boundary condition is obtained with $\gamma =\alpha \cdot 1j\cdot k_0$. For the excitation, an incident field $E_{Z,inc}=1$ is added with $q=\alpha \cdot 2j\cdot k_0$ on the left boundary. The natural (Neumann) boundary condition is a special case with $\gamma =q=0$, which does not have to be stated explicitly.

fem.assemble_boundaries_3rd(gamma={'bound_xmin': 1 / mu_plastic * 1j * k0, 
                                   'bound_xmax': 1 / mu_plastic * 1j * k0}, 
                            q={'bound_xmin': 1 / mu_plastic * 2j * k0, 
                               'bound_xmax': 0})

4. Solve the linear system for the field solution

fem.solve()

5. Process the solution

After solving, the field solution is stored in fem.phi as a complex vector. The individual real and imaginary parts are also stored seperately in fem.phi_re and fem.phi_im.

The corresponding locations of the N elements in the solution vector on the mesh are stored in fem.basis.doflocs, which has the shape (2, N).

scikit-fem offers helper functions to find certain elements, for example by means of labeled subdomains or boundaries:

x_bound_xmin, y_bound_xmin = fem.basis.doflocs[:, fem.basis.get_dofs('bound_xmin')]

Plotting is simple with the functions in skfem.visuals:

from skfem.visuals.matplotlib import plot
import matplotlib.pyplot as mplt

fig, ax = mplt.subplots(2, 1)
plot(fem.basis, fem.phi_re, ax=ax[0])
plot(fem.basis, fem.phi_im, ax=ax[1])
ax[0].set_aspect(1)
ax[1].set_aspect(1)
mplt.tight_layout()
mplt.show()