Mass matrix for a distributed vector problem

[bobsp...@gmail.com]

Dear all,

First of all, thanks for the awesome library. I am starting to learn deal.II, and it is great! The documentation has been extremely helpful so far.

I am trying to solve the time dependent elastic equation, and there we naturally have a mass matrix-like for the velocities. To build this mass matrix, I adapted the function Diffusion::assemble_system from step-52. However, in my problem the mass matrix has zero determinant, therefore, it is not invertible. Is there some mistake in the construction of the mass matrix, or something I am missing? 

I have sent attached a minimal source code that reproduces this error, the output of a single run, and a file containing the mass matrix.

I am running deal.II 9.0.0 installed via candi in Ubuntu 18.04.

Kind regards,

Bob

[luca.heltai]

You are integrating using two quadrature points per direction. Can you raise that to (2*fe.degree+1)?

L.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <matrix.txt><output.txt><minimal.cc>

[Robert Spartus]

Dear Luca,

Thanks for your suggestion. Unfortunately, it did not solve the problem. I am sending the modified version, as well as the output of the program.

Out of curiosity, what is the reason to use (2*fe_degree + 1)? Checking step-8, I notice that there a quadrature degree one larger than the polynomial degree is also used.

Bests,

Bob

[luca.heltai]

If you plan to use any domain that is not a square (or an affine transformation), you have to make sure you integrate exactly the product of two polynomials of order degree and of the determinant of the Jacobian. This last term is constant only for simple meshes, but it is the square root of a polynomial of order (degree-1) in more complicated cases.

2*fe_degree is ok for most cases, but I would not use this for serious calculations. I prefer to be on the safe side…

:)

L.

[Robert Spartus]

Dear Luca,

That is some fascinating information! It seems like step-44, for instance, does not follow this recommendation, as there the polynomial degree is 2, while the quadrature degree is 3, instead of the recommended 5 (https://dealii.org/developer/doxygen/deal.II/step_44.html#FiniteElementsystem). Is it because of the exception for squares you mentioned?

Do you have a reference that goes in depth on the topic of the choice of quadrature degrees? If so, I would grandly appreciate if you could send it my way.

Incidentally, have you been able to give any thought on the singularity of the mass matrix, even with the quadrature order is high?

Kind regards,

Bob

[Wolfgang Bangerth]

On 4/12/19 8:41 AM, Robert Spartus wrote:
>
> That is some fascinating information! It seems like step-44, for
> instance, does not follow this recommendation, as there the polynomial
> degree is 2, while the quadrature degree is 3

Actually, Gauss quadrature with degree+1 points in each direction is
sufficient to retain the convergence order of the finite element in
question, on any kind of mesh. Using higher order quadrature formulas
might increase the *absolute accuracy*, but is not necessary for the
convergence *order*.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[luca.heltai]

Wolfgang, is that true also for mass matrices? I’d agree with you for stiffness matrices, but I’d surprised this worked ok for mass matrices as well.

If so, I’ve always been over integrating in my life…

:)

L.

[Wolfgang Bangerth]

On 4/12/19 1:55 PM, luca.heltai wrote:
> Wolfgang, is that true also for mass matrices? I’d agree with you for
> stiffness matrices, but I’d surprised this worked ok for mass
> matrices as well.

I'm pretty sure. The theory goes like this: instead of computing the
matrix and rhs using the bilinear and linear forms

a(u,v) = f(v)

you're committing a variational crime by using quadrature instead of
integrals:

\tilde a(u,v) = \tilde f(v)

You then need to quantify the error due to this crime, and it turns out
that in order to not lose a convergence order, all you have to do is
compute the integrals via quadrature to the same convergence order as
for the overall finite element method. So, if you use elements of degree
k, you get O(h^k) in the energy norm, and you only need to integrate
matrix and rhs terms accurately enough to order O(h^k), which you can do
by using k+1 Gauss points in each coordinate direction.

[Robert Spartus]

Dear all,

Thanks for the insightful discussion on the integrating issue. Wolfgang, I guess your last argument is the same as you gave in one of your fantastic lectures?

Incidentally, do you have any ideas on how to solve the singularity of the mass matrix in this vector-valued problem?

Bests,

Bob

[Wolfgang Bangerth]

On 4/14/19 11:59 PM, Robert Spartus wrote:
>
> Thanks for the insightful discussion on the integrating issue. Wolfgang, I
> guess your last argument is the same as you gave in one of your fantastic
> lectures?

Yes. (Also, thanks for the compliment :-) )

> Incidentally, do you have any ideas on how to solve the singularity of the
> mass matrix in this vector-valued problem?

It is hard to imagine situations in which the mass matrix would be singular.
It is a positive definite form that gives rise to the mass matrix and so it
really shouldn't be singular at all. Can you show the code again with which
you build it?

The only situation where the mass matrix may become singular is if you have a
degenerate mesh with elements of zero or negative volume.

Best
WB

[Robert Spartus]

Dear Wolfgang,

> It is hard to imagine situations in which the mass matrix would be singular.
> It is a positive definite form that gives rise to the mass matrix and so it
> really shouldn't be singular at all. Can you show the code again with which
> you build it?

It seems that my mesh is neither singular nor degenerate. I wonder if the problem is that I am solving a vector valued problem. To build this mass matrix, I adapted the function Diffusion::assemble_system from step-52.

You will find the code attached, as well as the output of one run. If I isolate the mass matrix built and calculate its determinant in Python, the result is indeed zero.

I am really at a loss here, so any help is appreciated.

Bests,

Bob

[Robert Spartus]

Dear all,

I was doing some other investigations on this problem, and my findings made me even more confused.

I started analyzing the mass matrix generated by step-52. In this case, the only modification I did to the original was the following (full source attached): 

247,248d246

<     std::ofstream mat("matrix.txt");

<     mass_matrix.print_formatted(mat, 3, true, 0, "0.00", 1.0);

That is, we are writing the matrix into a file. All the rest of the program is unchanged.

Then, I used the following python code to calculate the determinant of this matrix (also attached):

     import numpy as np

     mat = np.loadtxt("matrix.txt")

     print("Determinant of mass matrix: {}".format(np.linalg.det(mat))) 

with output

    Determinant of mass matrix: 0.0.

This would mean that even in the *original* step-52 program the mass matrix is singular! However, this step runs without any issues. 

Does anyone have a piece of advice on how to solve this, and if I am making any mistakes?

Incidentally, you will find the matrix and a heat map of its non null entries attached.

Bests,

Bob

[Wolfgang Bangerth]

On 4/15/19 8:46 AM, Robert Spartus wrote:
>
> > It is hard to imagine situations in which the mass matrix would be singular.
> > It is a positive definite form that gives rise to the mass matrix and so it
> > really shouldn't be singular at all. Can you show the code again with which
> > you build it?
>
> It seems that my mesh is neither singular nor degenerate. I wonder if the
> problem is that I am solving a vector valued problem. To build this mass
> matrix, I adapted the function Diffusion::assemble_system from step-52.

I don't think the function you have gives you what you want. You have this:

cell_mass_matrix(i, j) +=
fe_values.shape_value(i, q_point) *
fe_values.shape_value(j, q_point) *
fe_values.JxW(q_point);

If you read the documentation of FEValues::shape_value(), you will see that
for vector-valued elements, it returns the one nonzero component of the vector
shape function. So shape_value(i)*shape_value(j) will always return something
nonzero. But what you really mean to do in your case is to multiply the
*vector* shape function i times the vector shape function j. Both of these
vectors may have a nonzero entry, but their dot product will only be nonzero
if these components are the same.

You will want to write the mass matrix here with extractors (i.e.,
fe_values[...]) in the same way you would build any other matrix for
vector-valued problems.

> You will find the code attached, as well as the output of one run. If I
> isolate the mass matrix built and calculate its determinant in Python, the
> result is indeed zero.

Using the determinant is an unreliable technique for large matrices. Think
about the case where you have a 1000x1000 matrix with eigenvalues all equal to
0.1. This is a perfectly invertible matrix, but the determinant is 0.1^1000,
which is zero for all practical purposes. You really need to look at the
eigenvalues themselves.

But, seeing the issue I mentioned above, if you have two components, then you
currently are computing the matrix
[ M M ]
[ M M ]
instead of
[ M 0 ]
[ 0 M ]
The former is clearly singular, so I'm not surprised.

Best
W.

[Robert Spartus]

Dear Wolfgang.

Thanks a ton for your input! You are completely right. Now that is you mention it, the problem is completely obvious. 

I really appreciate you pulling me out of that hole. 

Bests,

Bob

[Wolfgang Bangerth]

On 4/16/19 12:15 PM, Robert Spartus wrote:
>
> Thanks a ton for your input! You are completely right. Now that is you
> mention it, the problem is completely obvious.
>
> I really appreciate you pulling me out of that hole.

You're welcome. I'm glad it helped!