Quick update. I’ve moved to scipy
for easier matrix slicing and tried my luck at axing away the homogeneous dofs using this post to obtain total DoFs and taking its complementary to obtain a mask of free indices and that post to transform the PETSc Mat into a scipy CSR :
import dolfinx, ufl
import numpy as np
import scipy.sparse as sps
import scipy.sparse.linalg as la
from mpi4py.MPI import COMM_WORLD
n=100
mesh=dolfinx.UnitSquareMesh(COMM_WORLD,n,n)
FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V=dolfinx.FunctionSpace(mesh,FE)
u=ufl.TrialFunction(V)
v=ufl.TestFunction(V)
# Defining the forms
a=ufl.inner(ufl.grad(u),ufl.grad(v))*ufl.dx
b=ufl.inner( u, v )*ufl.dx
# Identify BC DoFs
dim = mesh.topology.dim - 1
mesh.topology.create_connectivity(dim, mesh.topology.dim)
boundary = np.where(np.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, dim, boundary)
# Total number of DoFs
N = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
# Free DoFs
freeinds = np.setdiff1d(range(N),boundary_dofs,assume_unique=True).astype(np.int32)
# Matrices
A = dolfinx.fem.assemble_matrix(a); A.assemble()
# Convert from PeTSc to Scipy for easy slicing
ai, aj, av = A.getValuesCSR()
# Slice away BCs
A = sps.csr_matrix((av, aj, ai))[freeinds,:][:,freeinds]
# Matrices
B = dolfinx.fem.assemble_matrix(b); B.assemble()
bi, bj, bv = B.getValuesCSR()
B = sps.csr_matrix((bv, bj, bi))[freeinds,:][:,freeinds]
# Solver
k=5
vals, vecs = la.eigs(A, k=k, M=B, sigma=0)
for i in range(k): print(np.abs(np.pi**2*((i//2)**2+(i//2+i%2)**2)+vals[i]))
I still get absurd errors, starting at 20
and going all the way up to 177
. Anyone has an idea out there ?