How to lump a mass matrix

1

Hi,

I am using Taylor-Hood element solving fluid problem with mesh motion (CG2 (velocity)-CG2(motion)-CG1(pressure)). To enhance stability, I want to do mass lumping. I am wondering if this is feasible in dolfin (Python).

Thanks in advance,
Victor

2

Victor,

Try the solution posted here: https://fenicsproject.org/qa/4284/mass-matrix-lumping/. I use lumping myself in elastodynamics.

Rudy.

3

Thanks Rudy,

I actually read this post a few hours ago, but noticed the last comment at the bottom on the page you gave saying “Be advised, however, that with parabolic and/or higher order continuous Lagrange interpolation lumping by using action (i.e. summing the elements of the consistent mass matrix rows) can lead lead to null or negative diagonal terms”.

So for your use case, are you using high order elements?

Thanks,
Victor

4

Victor,

I am not using higher order elements myself. According to Hughes’ “The Finite Element Method” textbook “no general theory of obtaining higher-order accurate mass matrices exists yet”. There is a wide range of special lumping techniques out there in the literature you may want to explore. In case you will not achieve a diagonalized mass matrix, the consistent mass is probably a good starting point.

Rudy.

5

Hi Victor,

did you find a way to implement the mass lumping for CG2 element? i’m also interested in that. Thanks.

Yang

6

Hi Yang,

Not really. But not needed anymore.

7

Hi,
you can look here:
https://comet-fenics.readthedocs.io/en/latest/demo/tips_and_tricks/mass_lumping.html

8

Hi Bleyer,

excellent! This is exactly what i’m looking for. Thanks a lot!

Additionally, I notice that for quadrilateral mesh, the singular matrix by “first” mass lumping way is no longer observed. Thus, i’m interested in comparing your implementation for quadrilateral mesh. Could you please provide some reference about how to set up the array “x” and “w” for quadrilateral mesh? Thank you!

Best regards,
Yang

9

Hello,
if you want to apply the same procedure for quadrilaterals, x correspond to the vertex positions in the reference element space and w are the corresponding weights for each points. Here I used 1D Gauss-Lobatto quadrature points and weights for 3 points:

elif isinstance(ref_el, UFCQuadrilateral)  and degree == 2:
    x = numpy.array([[0.0, 0.0],
               [0.5, 0.0],
               [1.0, 0.0],
               [0.0, 0.5],
               [0.5, 0.5],
               [1.0, 0.5],
               [0.0, 1.0],
               [0.5, 1.0],
               [1.0, 1.0]])
    w_1D = numpy.array([1/6., 2/3., 1/6.])
    w = numpy.kron(w_1D, w_1D)
    return QuadratureRule(ref_el, x, w)
1 Like
10

thank you so much, I will try it!