# Electromagnetic wave propagation in waveguide junction

I’m trying to solve the electromagnetic scattering problem with the initial wave TE10 coming from the left port of the waveguide junction. The weak form for the proble is shown in the following figure:

.
The problem has complex numbers and so I used the mixed function space to solve the problem. My code is:

from fenics import *
from mshr import *
from scipy.constants import c
import numpy as np

points_list = []
unit_coeff = 0.1
l = 11.17113576837
w = 11.43
initial_points = [[-l, w], [-l, -w], [-6.599113576837424, -w]]
initial_points = (np.array(initial_points)*unit_coeff).tolist()
[points_list.append(x) for x in initial_points]

for phi in [2 * np.pi/3, 4 * np.pi/3]:
rot_matrix = np.array([[np.cos(phi), -np.sin(phi)], [np.sin(phi), np.cos(phi)]])
for i, coord in enumerate(initial_points):
point = np.dot(rot_matrix, coord)
points_list.append(point.tolist())
points_list.append(points_list[0])
polygon_list = [Point(x, y) for x, y in points_list]

domain = Polygon(polygon_list)
circle = Circle(Point(0, 0), 0.35)
domain.set_subdomain(1, circle)
domain.set_subdomain(2, domain-circle)
mesh = generate_mesh(domain, 48)
materials = MeshFunction('size_t', mesh, 2, mesh.domains())

gm = 2.8e-3
f0 = gm * 200
fm = gm * 1317
f = 10
c = c * 100
fa = gm*135
c1 = 1 + fm * (f0+1j * fa) / ((f0+1j * fa) ** 2 - f**2)
c2 = -f*fm / ((f0+1j * fa) ** 2 - f**2)

class K1Real(UserExpression):
def __init__(self, materials, **kwargs):
super().__init__(degree=kwargs['degree'])
self.materials = materials
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = (c1 / (c1 ** 2 - c2 ** 2)).real
else:
values[0] = 1

class K1Imag(UserExpression):
def __init__(self, materials, **kwargs):
super().__init__(degree=kwargs['degree'])
self.materials = materials
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = (c1 / (c1 ** 2 - c2 ** 2)).imag
else:
values[0] = 0

class K2Real(UserExpression):
def __init__(self, materials, **kwargs):
super().__init__(degree=kwargs['degree'])
self.materials = materials
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = (c2 / (c1 ** 2 - c2 ** 2)).real
else:
values[0] = 0

class K2Imag(UserExpression):
def __init__(self, materials, **kwargs):
super().__init__(degree=kwargs['degree'])
self.materials = materials
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = (c2 / (c1 ** 2 - c2 ** 2)).imag
else:
values[0] = 0

class Epsilon(UserExpression):
def __init__(self, materials, **kwargs):
super().__init__(degree=kwargs['degree'])
self.materials = materials
def eval_cell(self, values, x, cell):
if self.materials[cell.index] ==1:
values[0] = 11.7
else:
values[0] = 1

g = Expression('cos(pi*x[1]/w)', degree=5, w=2 * w * 0.1, pi=np.pi)

tol = 0.0001
boundary_markers = MeshFunction('size_t', mesh, 1)

class PortOne(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], initial_points[0][0])

class PortTwo(SubDomain):
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.x0, self.y0 = points_list[3]
self.x1, self.y1 = points_list[4]
def inside(self, x, on_boundary):
return on_boundary and near((self.x1-self.x0) * (x[1]-self.y0) - (self.y1-self.y0) * (x[0]-self.x0), 0, tol)

class PortThree(SubDomain):
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.x0, self.y0 = points_list[6]
self.x1, self.y1 = points_list[7]
def inside(self, x, on_boundary):
return on_boundary and near((self.x1-self.x0) * (x[1]-self.y0) - (self.y1-self.y0) * (x[0]-self.x0), 0, tol)

class PEC(SubDomain):
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.terminal_points = [(points_list[1], points_list[2]),
(points_list[2], points_list[3]),
(points_list[4], points_list[5]),
(points_list[5], points_list[6]),
(points_list[7], points_list[8]),
(points_list[8], points_list[0])]
def inside(self, x, on_boundary):
on_pec_boundary = False
for pair_points in self.terminal_points:
x0, y0 = pair_points[0]
x1, y1 = pair_points[1]
on_pec_boundary = on_pec_boundary or near((x1-x0) * (x[1]-y0) - (y1-y0) * (x[0]-x0), 0, tol)
return on_boundary and on_pec_boundary

port_one = PortOne()
port_two = PortTwo()
port_three = PortThree()
pec_bc = PEC()

[b.mark(boundary_markers, i+1) for i, b in enumerate([port_one, port_two, port_three, pec_bc])]

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

P1 = FiniteElement('P', mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1 * P1)

u_r, u_i = TrialFunction(V)
v_r, v_i = TestFunction(V)

bc1 = DirichletBC(V.sub(0), Constant(0), boundary_markers, 4)
bc2 = DirichletBC(V.sub(1), Constant(0), boundary_markers, 4)
bcs = [bc1, bc2]

k1_r = K1Real(materials, degree=0)
k1_i = K1Imag(materials, degree=0)
k2_r = K2Real(materials, degree=0)
k2_i = K2Imag(materials, degree=0)
e_r = Epsilon(materials, degree=0)
k0 = 2 * np.pi * f / c * 1e9
gm = Constant((k0 ** 2 - (np.pi/(2 * w * 0.1))**2)**0.5)

class BaseComplexProduct(object):
def __init__(self):
self.r_r = 0
self.r_i = 0
def real(self):
return self.r_r + self.r_i
def imag(self):
return self.r_r - self.r_i

class ComplexDot(BaseComplexProduct):
def __init__(self, u_r, u_i, v_r, v_i):
self.r_r = dot( u_r, v_r ) + dot( u_i, v_i )
self.r_i = dot( u_i, v_r ) - dot( u_r, v_i )

class ComplexCrossProduct(BaseComplexProduct):
def __init__(self, ux_r, ux_i, uy_r, uy_i, vx_r, vx_i, vy_r, vy_i):
self.r_r = ux_r * vy_r + ux_i * vy_i - uy_r * vx_r - uy_i * vx_i
self.r_i = ux_i * vy_r - ux_r * vy_i - uy_i * vx_r + uy_r * vx_i

class ComplexProduct(BaseComplexProduct):
def __init__(self, u_r, u_i, v_r, v_i):
self.r_r = u_r * v_r + u_i * v_i
self.r_i = u_i * v_r - u_r * v_i

k2_r * ComplexCrossProduct(u_r.dx(0), u_i.dx(0), u_r.dx(1), u_i.dx(1),
v_r.dx(0),  v_i.dx(0), v_r.dx(1), v_i.dx(1)).imag() * dx + \
k2_i * ComplexCrossProduct(u_r.dx(0), u_i.dx(0), u_r.dx(1), u_i.dx(1),
v_r.dx(0),  v_i.dx(0), v_r.dx(1), v_i.dx(1)).real() * dx - \
e_r * pow(k0, 2) * ComplexProduct(u_r, u_i, v_r, v_i).real() * dx + \
gm * ComplexProduct(u_r, u_i, v_r, v_i).imag() * (ds(1) + ds(2) + ds(3))
L = 2 * gm* g * (v_r + v_i) * ds(1)

u = Function(V)
solve(a==L, u, bcs)

u_real, u_imag = u.split()

vtkfile = File('real.pvd')
vtkfile << u_real
vtkfile = File('imag.pvd')
vtkfile << u_imag


My question is that the solution is not symmetry about the x axis as indicated in the above figure, while the domain and the injection wave is symmetry about the x axis. Can anyone help me out. Thank you in advance!

My first guess is that this is what you get from a non-zero \kappa. Isn’t that how gyromagnetic elements like circulators work?
Just for me to understand: You’ve got an air-filled waveguide with a rod of \varepsilon_r = 11.7. \kappa_1 corresponds to the out-of-plane-component and \kappa_2 to the in-plane-component of the permeability tensor?

I tried to check your script by removing all non-essential stuff to get to a purely dielectrically loaded waveguide structure (removing all k1_r, k1_i, k2_r, k2_i in a, essentially setting \mu_r to 1 everywhere).
That makes for a more symmetric solution, which is what I would expect.

However, there I noticed that something may be off with your excitation g. It does not seem to come out as a proper cosine. It’s getting squished at the lower end. I did not quite understand why.

But there you are: Looks like you’ve built a working circulator with some oddities.