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.

Hello, how can I convert this code to Fenicsx/Dolfinx, for example Dolfinx 0.7.2 and Fenicsx 0.8, I wrote some codes but I faced with different errors. Also, if only want to obtain the solution of E field only for the fundemental mode (TE10) what should I do? Can you help to solve 2D equations for waveguides using Nedelec basis functions?

I converted in this way,

import numpy as np
import ufl
from petsc4py import PETSc
from slepc4py import SLEPc
from mpi4py import MPI
from dolfinx import mesh, fem
import dolfinx
import dolfinx.fem.petsc

Check for PETSc and SLEPc

if not hasattr(PETSc, ‘Mat’):
print(“PETSc is not available. Exiting.”)
exit()

if not hasattr(SLEPc, ‘EPS’):
print(“SLEPc is not available. Exiting.”)
exit()

Create mesh

width = 1.0
height = 0.5
mesh = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [width, height]], [8, 4], mesh.CellType.triangle)

Define function space

V = fem.FunctionSpace(mesh, (“Nedelec 1st kind H(curl)”, 2))

Define trial and test functions

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

Define the forms

s = ufl.inner(ufl.curl(v), ufl.curl(u)) * ufl.dx
t = ufl.inner(v, u) * ufl.dx
L = ufl.inner(ufl.as_vector((0.0, 0.0)), v) * ufl.dx

Define boundary condition

zero = PETSc.ScalarType((0, 0))
bc = fem.dirichletbc(zero, fem.locate_dofs_geometrical(V, lambda x: np.full(x.shape[1], True)), V)

Assemble the stiffness matrix (S) and mass matrix (T)

A = dolfinx.fem.petsc.assemble_matrix(fem.form(s), bcs=[bc])
A.assemble()
B = dolfinx.fem.petsc.assemble_matrix(fem.form(t), bcs=[bc])
B.assemble()

Convert assembled matrices to PETSc.Mat

A_petsc = A
B_petsc = B

Solve the eigensystem using SLEPc

eigensolver = SLEPc.EPS().create()
eigensolver.setOperators(A_petsc, B_petsc)
eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eigensolver.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)
eigensolver.solve()

Process and extract the eigenvalues

nconv = eigensolver.getConverged()
cutoff = None
count = 0

if nconv > 0:
for i in range(nconv):
lr, lc = eigensolver.getEigenvalue(i)
if lr > 1 and lc == 0:
cutoff = np.sqrt(lr)
_, _, rr, _ = eigensolver.getEigenpair(i)
u = fem.Function(V)
u.vector.setArray(rr)
u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
print(f"Eigenvalue {i}: {cutoff}")
count += 1
if count > 4:
break
else:
print(“No eigenvalues found.”)

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

Blockquote

however, I faced with errors, and can you check the code?


ValueError Traceback (most recent call last)
Cell In[33], line 34
32 s = ufl.inner(ufl.curl(v), ufl.curl(u)) * ufl.dx
33 t = ufl.inner(v, u) * ufl.dx
—> 34 L = ufl.inner(ufl.as_vector((0.0, 0.0)), v) * ufl.dx
36 # Define boundary condition
37 zero = PETSc.ScalarType((0, 0))

File /usr/local/lib/python3.10/dist-packages/ufl/measure.py:386, in Measure.rmul(self, integrand)
384 domain, = domains
385 elif len(domains) == 0:
→ 386 raise ValueError(“This integral is missing an integration domain.”)
387 else:
388 raise ValueError(“Multiple domains found, making the choice of integration domain ambiguous.”)

ValueError: This integral is missing an integration domain.

Thank you

Please do not duplicate posts. I am going to close this one in favor of Assembling Process.