Helmholtz eigenvalue problem in a rectangle

\[\begin{split} \mathbb{S}_u \begin{cases} \Omega = [0, L_x] \times [0, L_y] \\ u_{\text{D}}\vert_{\partial\Omega}=0 \\ f = 0 \\ \end{cases} \end{split}\]

import numpy as np
from dolfinx.fem import FunctionSpace
from lucifex.mesh import rectangle_mesh, mesh_boundary
from lucifex.solver import (
    evp, BoundaryConditions, OptionsSLEPc,
)
from lucifex.plt import plot_colormap, save_figure
from lucifex.pde.eigen import helmholtz

Lx = 1.0
Ly = 1.0
Nx = 32
Ny = 32
mesh = rectangle_mesh(Lx, Ly, Nx, Ny)
boundary = mesh_boundary(
        mesh, 
        {
            "left": lambda x: x[0],
            "right": lambda x: x[0] - Lx,
            "lower": lambda x: x[1],
            "upper": lambda x: x[1] - Ly,
        },
)
bcs = BoundaryConditions(('dirichlet', boundary.union, 0.0))
fs = FunctionSpace(mesh, ('P', 1))

nev = 10
slepc = OptionsSLEPc(
    eps_tol=1e-10,
    eps_target=0.0,
    eps_nev=nev,
    eps_ncv=50,
    eps_max_it=1000,
    eps_which='smallest_real',
)
u_solver = evp(helmholtz, bcs, slepc)(fs)
u_solver.solve()

sqrt = lambda x: np.sqrt(np.real(x)) if np.real(x) > 0 else None
k_values = [sqrt(-e) for e in u_solver.eigenvalues]
k_values

[None,
 None,
 None,
 4.444667393076721,
 7.034411761557017,
 7.034411761557114,
 8.900046596911208,
 9.967342410168534,
 9.967342410168753,
 11.361152904419761]

for k, u in zip(k_values, u_solver.eigenfunctions):
    if k is not None:
        fig, ax = plot_colormap(u, title=f'$u_{{k={k:.3f}}}$')
        save_figure(f'u(x,y)_k={k:.3f}', thumbnail=(k is k_values[-1]))(fig)