Rectangular waveguide problem chapter 34 Fenics book

![Screenshot_20200205-175120|250x500]

Hi guys,
My colleague wants me to run this code for him on my laptop, however the Fenics used in the book is very old syntax (pre 2013) and the code isn’t all there. It’s just snippets. Does anybody have the entire code for that problem in modern syntax or can translate it? I’m stuck I can’t do it. Thanks very much for help. It’s a rectangular waveguide problem ch34 first problem, and you analyse the cutoff frequency at the end I think (this isn’t my area, he just asked me to translate the code but I’m stuck…

I can’t see your screenshot, but do you mean this demo?

I’m very sorry about my screenshot not working! Unfortunately as a new user I can only post one picture. Here is the picture of the book page with the problem I’m looking at. The code is given in snippets with parts missing and I’m not finding it easy filling the blanks. I think the problem goes on for 2 or 3 pages. I’ll attach another picture in another reply in this thread to show you the problem I’m currently working on that is similar to the book problem so you better understand my predicament. Thank you very much for the quick reply.

Here is my problem that is similar to the book problem.

The demo I linked you to is the basic example of how to find the cutoff frequency in a rectangular waveguide. Hopefully you can figure out the syntax and where you need to go from there to apply it to your specific problem.

Thanks for the reply. Is there a way to change your boundary conditions in the code? I can’t seem to spot any.

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

This is brilliant. Thank you so much. I’ve showed it to my professor and he verified it. It works excellent thank you. Is there a way to do neumann conditions? I’m guessing I just have to change the DIRICHLETBC function?

Yes, just use natural assembly without a DirichletBC. Perhaps you should examine the FEniCS demos.

With regards to Maxwell’s equations for a beginner:

  • Take a lot of time to understand what’s going on when you solve a simple Poisson equation with various boundary conditions here.
  • Now take this knowledge and understand what is required to solve an eigenvalue problem.
  • Penultimately come to terms with the insidious problems with standard finite elements and the Maxwell operator. This is why we use the Nedelec element and not standard finite elements for the time harmonic Maxwell equations here.
  • Finally take lots of time to understand what this means for the boundary conditions you’re implementing, since they’re imposed in a tangential sense on the exterior boundary.
1 Like

Ah, so neumann bc are default in DOLFIN. So the only thing that changes is the mesh being used and a couple of values. Other than that I’d solve like I would for any other mesh?

There is no “default”. FEniCS’s philosophy is finite element matrix assembly of a specified weak formulation. I recommend the demos, and getting to grips with the basics of finite element theory.

Is there any way I could see an example of natural assembly for a rectangle waveguide so I can graph neumann conditions? I amn’t much of a programmer unfortunately.

Note that vtkplotter has been renamed to vedo, so just substitute

from vtkplotter.dolfin import plot

with

from vedo.dolfin import plot

Also, on Debian you may install it by running sudo apt install python3-vedo.