In this paper there’s a mixed formulation of Maxwells equations:
Find (\lambda, e, p)\in\mathbb{R}\times H_0(\operatorname{curl}; \Omega)\times H^1_0 (\Omega) such that e \neq d and:
(a) (\operatorname{curl} e, \operatorname{curl} \psi) + (\operatorname{grad}p, \psi)=\lambda(e, \psi),\quad\forall\psi\in H_0(\operatorname{curl}; \Omega)
(b) (e, \operatorname{grad} q)=0,\quad \forall q\in H_0^1(\Omega)
I’m trying to implement it in DolfinX like so
N1curl = ufl.FiniteElement("Nedelec 1st kind H(curl)", mesh.ufl_cell(), 1)
H1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = dolfinx.FunctionSpace(mesh, ufl.MixedElement(N1curl, H1))
e, p = ufl.TrialFunctions(V)
psi, q = ufl.TestFunctions(V)
a = ufl.inner(ufl.curl(e), ufl.curl(psi)) * ufl.dx + ufl.inner(ufl.grad(p), psi) * ufl.dx
a += ufl.inner(e, ufl.grad(q)) * ufl.dx
b = ufl.inner(e, psi) * ufl.dx
but I get a Zero pivot in LU factorization SLEPc error. What am I doing wrong?
Full code:
import dolfinx
import dolfinx.io
import ufl
from mpi4py import MPI
from slepc4py import SLEPc
import numpy as np
mesh = dolfinx.RectangleMesh(
MPI.COMM_WORLD,
[np.array([0, 0, 0]), np.array([np.pi, np.pi/2, 0])], [40, 40],
cell_type=dolfinx.cpp.mesh.CellType.triangle)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
N1curl = ufl.FiniteElement("Nedelec 1st kind H(curl)", mesh.ufl_cell(), 1)
H1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = dolfinx.FunctionSpace(mesh, ufl.MixedElement(N1curl, H1))
e, p = ufl.TrialFunctions(V)
psi, q = ufl.TestFunctions(V)
a = ufl.inner(ufl.curl(e), ufl.curl(psi)) * ufl.dx + ufl.inner(ufl.grad(p), psi) * ufl.dx
a += ufl.inner(e, ufl.grad(q)) * ufl.dx
b = ufl.inner(e, psi) * ufl.dx
bc_facets = np.where(np.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
bc_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, bc_facets)
u_bc = dolfinx.Function(V)
with u_bc.vector.localForm() as loc:
loc.set(0)
bc = dolfinx.DirichletBC(u_bc, bc_dofs)
A = dolfinx.fem.assemble_matrix(a, bcs=[bc])
A.assemble()
B = dolfinx.fem.assemble_matrix(b, bcs=[bc])
B.assemble()
eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
eps.setPowerShiftType(SLEPc.EPS.PowerShiftType.CONSTANT)
eps.setTarget(2)
eps.setDimensions(4)
eps.solve()
vals = [(i, eps.getEigenvalue(i)) for i in range(eps.getConverged())]
vals.sort(key=lambda x: x[1].real)
print(vals)
j = 0
E = dolfinx.Function(V)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "out.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
for i, _ in vals:
print(eps.getEigenpair(i, E.vector))
E.name = f"E-{j}"
xdmf.write_function(E)
j += 1