Problem with lumping, mixed function space and wrong lumped mass matrix!

Question

Dear all,
I have a very I have to solve a system of equation similiar to Cahn Hilliard model but I have two obstacles: first Mixed Function Space and second I must use lumped product scalar.
In particular, suppose that my weak form is :

(1/dt)*c*v*dx -(1/dt)* c0*v*dx + inner(grad(mu), grad(v)) = inner (grad( mu0), grad(v))
c*w*dx + relax*inner(grad(c), grad(w)) + {terms}  -relax* mu*w*dx = f*w*dx

where c, w are the variable to find, c0 the variable at the previous time step and f a known function. I MUST lump the terms

M11 =  (1/dt)*c*v*dx
M21 = c*w*dx
M22 = relax* mu*w*dx 

I use the method I've found in http://fenicsproject.org/qa/4284/mass-matrix-lumping
but I have problems with:
a21_lump = cwdx
mass_matrix21 = assemble(a21_lump)
mass_matrix21.zero()
mass_form21= action(a21_lump, const)
mass_matrix21.set_diagonal(assemble(mass_form21))
Indeed,it gives me (on a square mesh, 2x2)
[[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 1.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 4.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 4.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 3. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 9. 0. 0. 0. 0. 0. 0. 0. 0. ]
..... Other rows

INSTEAD OF:

[[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 1.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 4.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 4.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 3. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 9. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
.... Other rows

To obtain the last matrix (the right one ) I handle with the values of the matrix and thus I obtain memory errors if the mesh is bigger.
What can I do? Is there a way to obtain the right matrix as petsc?
I know that: the rows depends on the equation of the system (i.e the test function) and the fact that there is a value in a determined position depends on the fact that I'm dealing with c or mu. The trick is that, before the lumping, M21 is right.. it has the nonzero values in the right position!

Thanks a lot,
Cristina

Comments

Are c and w from the same FunctionSpace?

If not, it is not clear to me what mass lumping should mean.

If mass_matrix21 is a rectangular matrix, then setting the diagonal of that matrix does not make much sense either.

Maybe you should rethink your algorithm.

---

The algorithm was written by J. Barret: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.23.1474&rep=rep1&type=pdf

The matrix are square, I only copied only a part of them for saving space in my questions (there is written "etc" - et cetera :D )

And c and mu are defined as

V    = FunctionSpace(mesh,"Lagrange",1)
ME   = MixedFunctionSpace([V,V])
u= TrialFunction(ME)
c, mu = split(u)
{some code...}
u= Function(ME)
Time Loop

---

I see nothing wrong with lumped matrix. Entire $i-th$ row is lumped into the diagonal $(i, i)$ entry of the matrix as it should. How do you get the matrix that you presume is correct?

---

With:
ME = MixedFunctionSpace([V,V])
u = TrialFunction(ME)
c, mu = split(u)
v, w = TestFunctions(ME)

think about the lumped mass matrix for c and for mu..

M21= cwdx
M22 = muwdx

If I build the lumped matrix as in http://fenicsproject.org/qa/4284/mass-matrix-lumping,
I obtain the same matrix, that is wrong. I have to build manually the right matrix M21. I solved my problem building the petsc and then changing the values thanks to finalize_tensor = False