Lumping a mass matrix

Objectives

We discuss here how to lump a mass matrix.

Download sources

• Python script

• Jupyter notebook

Explicit dynamics simulations require the usage of lumped mass matrices i.e. diagonal mass matrices for which inversion can be done explicitly.

We show how to do this for Lagrange elements using the Gauss-Lobatto-Legendre quadrature rule, see also the Quadrature schemes tour. For the high order case in wave propagation, this setting defines the so-called spectral element method.

import numpy as np
from mpi4py import MPI
import ufl
import dolfinx.fem.petsc
from dolfinx import fem, mesh

domain = mesh.create_unit_square(
    MPI.COMM_WORLD, 1, 1, cell_type=mesh.CellType.quadrilateral
)

for order in range(1, 3):
    V = fem.functionspace(domain, ("P", order))
    v = ufl.TestFunction(V)
    u = ufl.TrialFunction(V)

    dx = ufl.Measure("dx", domain=domain)
    dx_lumped = dx(metadata={"quadrature_rule": "GLL", "quadrature_degree": order})
    mass_form = v * u * dx
    lumped_mass_form = v * u * dx_lumped

    M_consistent = fem.assemble_matrix(fem.form(mass_form))
    print(
        "Consistent mass matrix:\n", np.array_str(M_consistent.to_dense(), precision=3)
    )

    M_lumped = fem.assemble_matrix(fem.form(lumped_mass_form))
    print("Lumped mass matrix:\n", np.array_str(M_lumped.to_dense(), precision=3))

INFO: Disabling color, you really want to install colorlog.

Consistent mass matrix:
 [[0.111 0.056 0.056 0.028]
 [0.056 0.111 0.028 0.056]
 [0.056 0.028 0.111 0.056]
 [0.028 0.056 0.056 0.111]]

Lumped mass matrix:
 [[0.25 0.   0.   0.  ]
 [0.   0.25 0.   0.  ]
 [0.   0.   0.25 0.  ]
 [0.   0.   0.   0.25]]

Consistent mass matrix:
 [[ 0.018 -0.004 -0.004  0.001  0.009  0.009 -0.002 -0.002  0.004]
 [-0.004  0.018  0.001 -0.004  0.009 -0.002  0.009 -0.002  0.004]
 [-0.004  0.001  0.018 -0.004 -0.002  0.009 -0.002  0.009  0.004]
 [ 0.001 -0.004 -0.004  0.018 -0.002 -0.002  0.009  0.009  0.004]
 [ 0.009  0.009 -0.002 -0.002  0.071  0.004  0.004 -0.018  0.036]
 [ 0.009 -0.002  0.009 -0.002  0.004  0.071 -0.018  0.004  0.036]
 [-0.002  0.009 -0.002  0.009  0.004 -0.018  0.071  0.004  0.036]
 [-0.002 -0.002  0.009  0.009 -0.018  0.004  0.004  0.071  0.036]
 [ 0.004  0.004  0.004  0.004  0.036  0.036  0.036  0.036  0.284]]

Lumped mass matrix:
 [[0.028 0.    0.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.028 0.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.028 0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.028 0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.111 0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.111 0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.111 0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.111 0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.    0.444]]

For explicit dynamics simulation, the mass matrix can then be manipulated using the diagonal vector, to compute for instance \(M^{-1}w\) where \(w\) is some Function.

mass_diagonal = u * dx_lumped  # defines a linear form corresponding to the diagonal
M_vect = fem.petsc.assemble_vector(fem.form(mass_diagonal))
w = fem.Function(V)
iMw = fem.Function(V)
iMw.vector.pointwiseDivide(M_vect, w.vector)