I am trying to reproduce this tutorial about stable and unstable finite elements for Maxwell’s equations in Fenicsx. However, I am struggling with the set up of the SLEPc solver. So far, I have the following MWE, where only the Lagrangian elements are working yet:

```
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from slepc4py import SLEPc
from dolfinx import mesh
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, form, dirichletbc, locate_dofs_geometrical, locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix
from ufl import FiniteElement, VectorElement, TestFunction, TrialFunction, curl, dot, dx, grad, inner
def eigenvalues(V, bcs):
u = TrialFunction(V)
v = TestFunction(V)
a = form(inner(curl(u), curl(v))*dx)
b = form(inner(u, v)*dx)
A = assemble_matrix(a, bcs)
A.assemble()
B = assemble_matrix(b, bcs)
B.assemble()
# [bc.zero(B) for bc in bcs] # in the fenics tutorial they zero out the corresponding lines, why?
solver = SLEPc.EPS()
solver.create()
solver.setOperators(A, B)
solver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
solver.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
solver.setTarget(5.5)
neigs = 12
solver.setDimensions(nev=neigs, ncv=2*neigs)
solver.setFromOptions()
solver.solve()
print(solver.getConverged())
# Return the computed eigenvalues in a sorted array
computed_eigenvalues = []
for i in range(min(neigs, solver.getConverged())):
r = solver.getEigenvalue(i) # ignore the imaginary part
computed_eigenvalues.append(r)
return np.sort(np.array(computed_eigenvalues))
ncells = 40
corners = [[0.0, 0.0], [np.pi, np.pi]]
msh = mesh.create_rectangle(MPI.COMM_WORLD, corners, [ncells, ncells])
V = VectorFunctionSpace(msh, ("CG", 1))
# V = FunctionSpace(msh, ("N1curl", 1))
def top(x):
return np.logical_and(
np.isclose(x[1], np.pi), x[0] < np.pi
)
def bottom(x):
return np.logical_and(
np.isclose(x[1], 0.0), x[0] > 0
)
def right(x):
return np.logical_and(
np.isclose(x[0], np.pi), x[1] < np.pi
)
def left(x):
return np.logical_and(
np.isclose(x[0], 0.0), x[1] > 0
)
bd_markers_y = [top, bottom]
bd_markers_x = [left, right]
bcs = []
for bd_marker in bd_markers_y:
bd_facets = mesh.locate_entities_boundary(msh, msh.topology.dim-1, bd_marker)
bd_dofs_x = locate_dofs_topological(V.sub(0), msh.topology.dim-1, bd_facets)
bc = dirichletbc(ScalarType(0), bd_dofs_x, V.sub(0))
bcs.append(bc)
for bd_marker in bd_markers_x:
bd_facets = mesh.locate_entities_boundary(msh, msh.topology.dim - 1, bd_marker)
bd_dofs_y = locate_dofs_topological(V.sub(1), msh.topology.dim - 1, bd_facets)
bc = dirichletbc(ScalarType(0), bd_dofs_y, V.sub(1))
bcs.append(bc)
ev = eigenvalues(V, bcs)
print(ev)
```

More precisely, I have two questions

- In the tutorial, the lines corresponding to the BCs in the operator B are zeroed out. How can this be done in Fenicsx?
- When I comment out the lines

```
solver.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
solver.setTarget(5.5)
```

I get the 12 eigenvalues with the largest magnitude as expected. However, when I include the lines, no eigenvalue converges. What could be the cause of this?