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
a = k1_r * ComplexDot(grad(u_r), grad(u_i), grad(v_r), grad(v_i)).real() * dx + \
k1_i * ComplexDot(grad(u_r), grad(u_i), grad(v_r), grad(v_i)).imag() * dx - \
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!
