Rectangular waveguide problem chapter 34 Fenics book

It’s been a long time since I’ve looked at waveguide problems. But here’s an example with assumed perfect conductor BCs

You need vtkplotter if you want the plotting, installed via pip3 install vtkplotter --user. Otherwise just remove the plotting lines. Just plotting the first 4 propagating modes.

Also I’m not sure how correct this is. You should construct a test with a solution known a priori to make sure you’re solving what you think you’re solving and don’t use this example as a black box.

from dolfin import *

width = 1.0
height = 0.5
mesh = RectangleMesh(Point(0, 0), Point(width, height), 8, 4, "left/right")

V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 2)

v = TestFunction(V)
u = TrialFunction(V)

s = inner(curl(v), curl(u))*dx
t = inner(v, u)*dx
L = dot(Constant((0.0, 0.0)), v)*dx

# Hermitian assembly
perf_conductor_bc = DirichletBC(V, Constant((0, 0)), "on_boundary")
S = PETScMatrix()
T = PETScMatrix()
dummy = PETScVector()
assemble_system(s, L, A_tensor=S, b_tensor=dummy, bcs=perf_conductor_bc)
assemble_system(t, L, A_tensor=T, b_tensor=dummy, bcs=perf_conductor_bc)

esolver = SLEPcEigenSolver(S, T)
esolver.parameters["spectrum"] = "smallest real"
esolver.parameters["solver"] = "lapack"
esolver.solve()

from vtkplotter.dolfin import plot
cutoff = None
count = 0
for i in range(S.size(1)):
    (lr, lc, rr, rc) = esolver.get_eigenpair(i)
    u = Function(V)
    if lr > 1 and lc == 0:
        cutoff = sqrt(lr)
        u.vector()[:] = rr

        plot(u, shape=((2,2)), at=count)
        count += 1
        if count > 4:
            break

if cutoff is None:
    print("Unable to find dominant mode")
else:
    print("Cutoff frequency:", cutoff)

2 Likes