Solving generalized eigenvalues of mixed finite element stiff matrix

Question

Hi friends,
I am considering to solve generalized eigenvalues of stiff matrix formed by mixed finite element, which involves five subspaces. But an error "Zero pivot row 9 value 0 tolerance 2.22045e-14!" puzzled me. I tried a lot but failed to resolve the error. My code is as follows:

Load mesh and subdomains

mesh = UnitSquareMesh(4,4)

Define function spaces

Vy = TensorFunctionSpace(mesh, "DG", 1, symmetry=True)
Vr = FunctionSpace(mesh, "CG", 2)
Vp = FunctionSpace(mesh, "CG", 1)
Vu = FunctionSpace(mesh, "CG", 2)
Vqp = VectorFunctionSpace(mesh, "CG", 2)
W = MixedFunctionSpace([Vy, Vr, Vp, Vu, Vqp])

Define boundary conditions

u0 = Constant(0.0)
u1 = Constant((0.0,0.0))
bc1 = DirichletBC(W.sub(1), u0, "on_boundary")
bc3 = DirichletBC(W.sub(3), u0, "on_boundary")
bc4 = DirichletBC(W.sub(4), u1, "on_boundary")
bcs = [bc1, bc3, bc4]

Define variational problem

(y, r, p, u, p1) = TrialFunctions(W)
(z, s, q, v, eta) = TestFunctions(W)

a = (inner(y, z) - inner(nabla_grad(p1),z) - dot(nabla_grad(s), nabla_grad(u)) \
+dot(nabla_grad(s), p1) + qrot(p1)-dot(nabla_grad(r), nabla_grad(v)) \
-inner(y, nabla_grad(eta))+dot(nabla_grad(r), eta) + prot(eta))*dx

m = uvdx

Assemble stiffness matrix

A = PETScMatrix()
assemble(a, tensor=A)
M = PETScMatrix()
assemble(m, tensor=M)

Apply the boundary conditions

for bc in bcs:
bc.apply(A)

bc3.apply(M)

Create eigensolver

eigensolver = SLEPcEigenSolver(A,M)
eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.parameters['tolerance'] = 1e-6

Solve for eigenvalues

eigensolver.solve()

w = Function(W)
for i in range(eigensolver.get_number_converged()):
#extract next eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(i)
print ("eigenvalue:", i, r, c)
#assign eigenvector to function
w.vector()[:] = rx
y, r, p, u, p1 = w.split()
#plot eigenfunction
plot(u)
interactive()

Does anybody can give me a hint? Thanks in advance!

Comments

Are you sure, that your form a is correct?

---

Yes, the weak form is correct. The problem is there are zero rows in M.

Answer 1

The matrix $M$ on the r.h.s. of the generalized eigenproblem is singular. In fact it is non-zeros only for the $u$-block, and zeros for the $y$, $r$, $p$, $p1$-blocks.

I believe that SLEPcEigenSolver only supports generalized eigenproblems with non-singular r.h.s. matrix.

Best,

Umbe

Comments

Umber, you are right. I output the matrices A and M, Matlab can work for me. Thanks.