# Mixed formulation of Maxwell's equations

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:

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


I know this doesn’t answer your question, but perhaps it’s a useful first step:

By the nature of the generalised eigenvalue problem you shouldn’t need to implement the involution term as you would for a time harmonic Maxwell problem where the frequency is not an eigenvalue of the system.

Consider the following. Let (u, \lambda > 0) solve the generealised eigenvalue problem, then we have

\begin{align} \nabla \times \nabla \times u &= \lambda u, \\ \nabla \cdot (\nabla \times \nabla \times u) &= \lambda \nabla \cdot u, \\ 0 &= \nabla \cdot u. \end{align}

If you modify your problem do you have a full diagonal and as such does the eigensolver converge?

Also there’s a very old (but excellent) demo here which may help.

Yes, I’ve seen this example. For some reason I don’t get the same eigenvalues as in the example.

According to this paper, presented in the question formulation, should better impose divergence free condition \nabla\cdot u = 0.

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, 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)
V = dolfinx.FunctionSpace(mesh, N1curl)

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

a = ufl.inner(ufl.curl(u), ufl.curl(v)) * ufl.dx + ufl.inner(ufl.curl(u), ufl.curl(v)) * ufl.dx
b = ufl.inner(u, v) * 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(5.5)
eps.setDimensions(12)
eps.solve()

np.isclose(eps.getEigenvalue(i), 1)]
vals = [(i, eps.getEigenvalue(i)) for i in range(eps.getConverged())]

vals.sort(key=lambda x: x[1].real)

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


For some reason I don’t get the same eigenvalues as in the example.

Wops, there’s a typo a = ufl.inner(ufl.curl(u), ufl.curl(v)) * ufl.dx. Example works fine.

Turns out, that the problem was not with the formulation, but with the configuration of the Spectral Transformation (ST) if SLEPc. It shold be configured as follows:

st = eps.getST()
st.setType(SLEPc.ST.Type.SINVERT)


Checking \int_{\Omega}\psi\;\text{d}x values I can see that unphysical modes with \lambda=1 all have high integral value:

eigenvalue: 0.999689889049;     int(psi): -6.308888402483e-14
eigenvalue: 0.999967476453;     int(psi): -1.645503611179e-12
eigenvalue: 1.000000000000;     int(psi): 1.337930720732e-03
eigenvalue: 1.000000000000;     int(psi): -1.598505351224e-03
eigenvalue: 1.000000000000;     int(psi): 2.309336061885e-02
eigenvalue: 1.000000000000;     int(psi): -1.404544335133e-03
eigenvalue: 1.000000000000;     int(psi): 3.885881138543e-04
eigenvalue: 1.000000000000;     int(psi): 6.228103117720e-04
eigenvalue: 1.000000000000;     int(psi): 1.856027911829e-03
eigenvalue: 1.000000000000;     int(psi): 3.188035050973e-03
eigenvalue: 2.000342166381;     int(psi): 4.142071959184e-14
eigenvalue: 3.997258892126;     int(psi): -4.388551693095e-14
eigenvalue: 3.997260387796;     int(psi): -1.283612474464e-14
eigenvalue: 4.997207026787;     int(psi): 7.494442252868e-15
eigenvalue: 5.002446610361;     int(psi): 1.745873985960e-14
eigenvalue: 8.005430745727;     int(psi): -2.022234716329e-14
eigenvalue: 8.984888327100;     int(psi): 4.700200964440e-14
eigenvalue: 8.987372947191;     int(psi): -5.860825488704e-14
eigenvalue: 9.992103624262;     int(psi): -2.649171526890e-14
eigenvalue: 9.992163510773;     int(psi): 1.579640268136e-14


Full DolfinX 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, 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)
st = eps.getST()
st.setType(SLEPc.ST.Type.SINVERT)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
eps.setTarget(5.5)
eps.setDimensions(20)
eps.solve()

vals = [(i, eps.getEigenvalue(i)) for i in range(eps.getConverged())]
vals.sort(key=lambda x: x[1].real)

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:
E1, E2 = E.split()
print(f"eigenvalue: {eps.getEigenpair(i, E.vector).real:.12f};\t"
f"int(psi): {dolfinx.fem.assemble_scalar(E2 * ufl.dx):.12e}")
E1.name = f"E1-{j:03d}-{eps.getEigenpair(i, E.vector).real:.4f}"
E2.name = f"E2-{j:03d}-{eps.getEigenpair(i, E.vector).real:.4f}"
xdmf.write_function(E1)
xdmf.write_function(E2)
j += 1

1 Like

I am doing something similar and I was wondering what the explanation for the unphysical modes with \lambda = 1 was?

1 Like

If you zero out the rows corresponding to the boundary conditions in the operator B, you actually do not retrieve spurious eigenmodes with eigenvalue 1.0, i.e.

B = assemble_matrix(b, bcs=[bc], diagonal=0.0)

1 Like

Hi, when I want to run this program, I get the error "AttributeError: module ‘dolfinx’ has no attribute ‘RectangleMesh’ "… Do you have any idea why? Thanks…

The code here is written for an older version of DOLFINx.
dolfinx.RectangleMesh has been replaced with dolfinx.mesh.create_rectangle.

Similarly, You should replace

with

dolfinx.fem.dirichletbc


and

with

A = dolfinx.fem.petsc.assemble_matrix(a, bcs=[bc])
A.assemble()
B = dolfinx.fem.petsc.assemble_matrix(b, bcs=[bc])
B.assemble()


Thank you very much for your help. In the meantime, I tried to reinstall everything related Dolfinx. Then I followed your fixes. Now I get the error " ModuleNotFoundError: No module named ‘basix.ufl_wrapper’ ". Searched for it but couldn’t find anything. Could you possibly help please?

This probably has to do with having non-compatible versions of Basix, ffcx and dolfinx. How did you install each individual component?

I tried many things. spack installation, ubuntu ppa, and individually (e.g., pip for ffcx). My system is Ubuntu 22.04.2 LTS. Can you suggest the best way to install everything as simple as it can be?

The problem is if you have mixed all of these

There can be many conflicting versions/dependencies from the different installation methods.

The easiest is to use conda.
Second easiest is to use docker
Third easiest is spack
Fourth easiest is ubuntu PPA(But not recommended as the package there is not up to date)

However, i cannot help you any further if you dont post a minimal code example that yields an error, and the full trace of the error message.

Many thanks. I will try reinstallation using conda. If the next try fails too, I will post the code with the error message.

If you decide to build everything manually (e.g. for the computational cluster), here’s a semi-automatic script

SOURCE_DIR=$(realpath .) BUILD_DIR=$SOURCE_DIR/build-dir
INSTALL_DIR=/opt/dolfinx

pip3 uninstall -y --break-system-packages fenics-dolfinx fenics-ffcx fenics-ufl fenics-basix

sudo apt install libpetsc-complex-dev libslepc-complex-dev python3-petsc4py-complex python3-slepc4py-complex

BASIX_COMMIT=185dd46f4b0121539fbc126ed338836f453b8c17
DOLFIN_COMMIT=3784623b46081f408052b42123bb290da158a3e9
UFL_COMMIT=6eeef938c877634dcb7c8acb101bb2363ed1455b
FFCX_COMMIT=c1dc25e7e6ec2661743449a27906f203760c0627

git clone --depth 1 https://github.com/KarypisLab/GKlib.git
git clone --depth 1 https://github.com/KarypisLab/METIS.git
git clone --depth 1 https://github.com/KarypisLab/ParMETIS.git $SOURCE_DIR/ParMETIS git clone --depth 1 https://github.com/ornladios/ADIOS2.git$SOURCE_DIR/ADIOS2
git clone https://github.com/FEniCS/basix.git $SOURCE_DIR/basix git -C$SOURCE_DIR/basix checkout $BASIX_COMMIT git clone https://github.com/FEniCS/dolfinx.git$SOURCE_DIR/dolfinx
git -C $SOURCE_DIR/dolfinx checkout$DOLFINX_COMMIT
git clone https://github.com/FEniCS/ufl.git $SOURCE_DIR/ufl git -C$SOURCE_DIR/ufl checkout $UFL_COMMIT git clone https://github.com/FEniCS/ffcx.git$SOURCE_DIR/ffcx
git -C $SOURCE_DIR/ffcx checkout$FFCX_COMMIT

scalar_type=complex

CMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH:$INSTALL_PATH/lib/cmake/:$INSTALL_DIR/share/cmake/" PETSC_DIR="/usr/lib/petscdir/petsc-$scalar_type"
SLEPC_DIR="/usr/lib/slepcdir/slepc-$scalar_type" PETSC_ARCH="linux-gnu-$scalar_type-64"
PYTHONPATH="$PETSC_DIR/lib/python3/dist-packages:$SLEPC_DIR/lib/python3/dist-packages"
PKG_CONFIG_PATH="$PETSC_DIR/lib/pkgconfig:$PKG_CONFIG_PATH"
LD_LIBRARY_PATH="$PETSC_DIR/lib:$LD_LIBRARY_PATH"

# GKlib
pushd GKlib
make config prefix=$INSTALL_DIR sudo make install popd # METIS pushd METIS make config prefix=$INSTALL_DIR
sudo make install
popd

# ParMETIS
pushd ParMETIS
make config cc=mpicc prefix=$INSTALL_DIR sudo make install popd # ADIOS2 rm -rf$BUILD_DIR
mkdir -p $BUILD_DIR cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALL_DIR \
-B $BUILD_DIR -S$SOURCE_DIR/ADIOS2
cmake --build $BUILD_DIR -j$(nproc)
sudo cmake --install $BUILD_DIR # Basix rm -rf$BUILD_DIR
mkdir -p $BUILD_DIR cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALL_DIR \
-B $BUILD_DIR -S$SOURCE_DIR/basix/cpp
cmake --build $BUILD_DIR -j$(nproc)
sudo cmake --install $BUILD_DIR pushd$SOURCE_DIR/basix
pip3 install --break-system-packages --user .
popd

# UFL
pushd $SOURCE_DIR/ufl pip3 install --break-system-packages --user . popd # FFCX pushd$SOURCE_DIR/ffcx
pip3 install --break-system-packages --user .
popd

# DolfinX
rm -rf $BUILD_DIR mkdir -p$BUILD_DIR

cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALL_DIR \ -B$BUILD_DIR -S $SOURCE_DIR/dolfinx/cpp cmake --build$BUILD_DIR -j$(nproc) sudo cmake --install$BUILD_DIR
source $INSTALL_DIR/lib/dolfinx/dolfinx.conf pushd$SOURCE_DIR/dolfinx/python
pip3 install --break-system-packages --user .
popd


Note, that the same project tags are not necessarily compatible, that’s why I use commit hashes.

Thank you very much. I will try this as well.

Sorry again for disturbing. The above code with your fixes (given below) now results in the error:

in <module>
mesh = dolfinx.mesh.create_rectangle(
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/mesh.py", line 372, in create_rectangle
mesh = _cpp.mesh.create_rectangle(comm, points, n, cell_type, ghost_mode, partitioner, diagonal)
TypeError: create_rectangle(): incompatible function arguments. The following argument types are supported:
1. (comm: MPICommWrapper, p: List[List[float[2]][2]], n: List[int[2]], celltype: dolfinx.cpp.mesh.CellType, ghost_mode: dolfinx.cpp.mesh.GhostMode, partitioner: Callable[[MPICommWrapper, int, int, dolfinx::graph::AdjacencyList<long>, dolfinx.cpp.mesh.GhostMode], dolfinx::graph::AdjacencyList<int>], diagonal: dolfinx.cpp.mesh.DiagonalType) -> dolfinx.cpp.mesh.Mesh

Invoked with: <mpi4py.MPI.Intracomm object at 0x7ff8947af1e0>, [array([0, 0, 0]), array([3.14159265, 3.14159265, 0.        ])], [40, 40], <CellType.triangle: 3>, <GhostMode.shared_facet: 1>, <built-in method  of PyCapsule object at 0x7ff89232c1e0>, <DiagonalType.right: 1>



The code is exactly this:

import dolfinx
import dolfinx.io
import ufl
from mpi4py import MPI
from slepc4py import SLEPc
import numpy as np

mesh = dolfinx.mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([0, 0, 0]), np.array([np.pi, np.pi, 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.fem.dirichletbc(u_bc, bc_dofs)

A = dolfinx.fem.petsc.assemble_matrix(a, bcs=[bc])
A.assemble()
B = dolfinx.fem.petsc.assemble_matrix(b, bcs=[bc])
B.assemble()

eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
st = eps.getST()
st.setType(SLEPc.ST.Type.SINVERT)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
eps.setTarget(5.5)
eps.setDimensions(20)
eps.solve()

vals = [(i, eps.getEigenvalue(i)) for i in range(eps.getConverged())]
vals.sort(key=lambda x: x[1].real)

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:
E1, E2 = E.split()
print(f"eigenvalue: {eps.getEigenpair(i, E.vector).real:.12f};\t"
f"int(psi): {dolfinx.fem.assemble_scalar(E2 * ufl.dx):.12e}")
E1.name = f"E1-{j:03d}-{eps.getEigenpair(i, E.vector).real:.4f}"
E2.name = f"E2-{j:03d}-{eps.getEigenpair(i, E.vector).real:.4f}"
xdmf.write_function(E1)
xdmf.write_function(E2)
j += 1


Change this to
array([0, 0]), array([3.14159265, 3.14159265 ])

The error now is:
V = dolfinx.FunctionSpace(mesh, ufl.MixedElement(N1curl, H1))
AttributeError: module ‘dolfinx’ has no attribute ‘FunctionSpace’

dolfinx.fem.FunctionSpace.
I would strongly suggest that you familiarize yourself with the new api either by looking at the demos
https://docs.fenicsproject.org/dolfinx/v0.6.0/python/demos.html
or the tutorial
https://jsdokken.com/dolfinx-tutorial/