MultiMesh in solving eigenvalue problem of Stokes operator

[xfliu]

I am following your code to use multimesh to solve eigenvalue problem of the Stokes operator.

By selecting test function from a rawer mesh, I want to weaken the divergence-free condition.
When mesh_u and mesh_div_free are the same mesh, the eigenfunction should be completely divergence-free.
But the code below failed to apply the divergence-free condition to solution u.

Could you give me any hint where the code is wrong?

from fenics import *

mesh_u=UnitSquareMesh(32, 32, diagonal="right")
mesh_div_free=UnitSquareMesh(16, 16, diagonal="right")
mesh_R=UnitSquareMesh(1, 1, diagonal="right")

multimesh = MultiMesh()
multimesh.add(mesh_u)
multimesh.add(mesh_div_free)
multimesh.add(mesh_R)
multimesh.build()

degree=3

CG = VectorElement("CG", triangle, degree,2)
DG = FiniteElement("DG", triangle, degree-1)
Real = FiniteElement("Real", triangle, 0)
MixedV = MultiMeshFunctionSpace(multimesh, MixedElement([CG, DG, Real]))
V = MultiMeshSubSpace(MixedV, 0)

(u, p, alpha) = TrialFunctions(MixedV)
(v, q, beta) = TestFunctions(MixedV)

a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx + alpha*q*dx + beta*p*dx
b = inner(u,v)*dx

f = Constant((0.0, 0.0))
L = inner(f,v)*dx 


mvc = MeshValueCollection("size_t", mesh_u, 1)
mf = cpp.mesh.MeshFunctionSizet(mesh_u, mvc)

noslip = Constant((0.0, 0.0))
marker=100
bc0=MultiMeshDirichletBC(V, noslip, mf, marker, 0)

A = assemble_multimesh(a)
B = assemble_multimesh(b)

bc0.apply(A)
bc0.zero(B)

# downcast to PETSc matrices
MatA = as_backend_type(A)
MatB = as_backend_type(B)

import scipy.io as myio
import scipy.sparse as sp

row_, col_, val_ = MatA.mat().getValuesCSR()
mat_A = sp.csr_matrix((val_,col_,row_))
row_,col_,val_ = MatB.mat().getValuesCSR()
mat_B = sp.csr_matrix((val_,col_,row_))

#Another approach to remove boundary value DOF
bd_int = inner(u,v)*ds
bd_int = assemble_multimesh(bd_int)
Mat_bd_int = as_backend_type(bd_int)

row_,col_,val_ = Mat_bd_int.mat().getValuesCSR()
Mat_bd_int = sp.csr_matrix((val_,col_,row_))
myio.savemat("./matrix/matrix.mat",{'A':mat_A,'B':mat_B,'BD_Int':Mat_bd_int})

(I exported the matrix to mat file and then solve it in Octave or MATLAB.)