Mass Lumping in case of Dirichlet boundary conditions

Question

I'm currently trying to implement a FEM to solve a type of wave equation with homogeneous Dirichlet boundary conditions by using standard $\mathcal{P}_1$ triangle elements and an explicit scheme for the time discretization. In every discrete time step, the solution of a linear system with the mass matrix $$M_{i,j} = \int_{\Omega}\varphi_i(x) \varphi_j(x) \, dx$$ occurs, where $i,j \in \mathcal{N_{\mathrm{inner}}}$ are the inner nodes and $\varphi_i$ are the piecewise linear, global continuous, nodal basis functions. Thus, i would like to replace $M$ by a good diagonal approximation $L$, namely the $\textit{lumped mass matrix}$. Now, i've encountered two types to lumping techniques. First, a row-sum technique $$L_{i,j} = \delta_{i, j} \sum_{k \, \in \, \mathcal{N_{\mathrm{inner}}}} M_{i,k}.$$ And second, a nodal quadrature technique $$L_{i,j} = \int_{\Omega} I_h(\varphi_i(x) \varphi_j(x)) \, dx,$$ where $I_h$ is the piecewise linear interpolator. In case of $\mathcal{P_1}$-elements, the above can be rewritten as $$L_{i,j} = \delta_{i, j} \sum_{T \, \in \, \omega(i)} \frac{|T|}{3},$$ with $\omega(i)$ being the patch of elements that share the inner node $i \in \mathcal{N_{\mathrm{inner}}}$.

In case of homogeneous Neumann boundary conditions, we replace $\mathcal{N_{\mathrm{inner}}}$ by the set $\mathcal{N}$ of all nodes on the mesh. Then, both approaches are the same, because $\sum_{i \, \in \, \mathcal{N}} \varphi_i(x) \equiv 1$ on $\bar{\Omega}$. But for homogeneous Dirichlet boundary conditions, this is not the case anymore and both techniques result in different matrices (though, the two lumped diagonal matrices are only different on nodes $i \in \mathcal{N_{\mathrm{inner}}}$ that are nodes of an element, which touches the boundary).

$\textbf{Now my Question}$: which one of these lumped mass matrices should I use? I've seen in some applications, that people want—for conservation reasons—that the row sum of both the original mass matrix and the lumped mass matrix are equal, but i've also seen that people want the node quadrature approach, because, in least in theory, the resulting fully discrete scheme shouldn't lose any order of accuracy.

Answer 1

Since you are using simple elements and an explicit solver, my approach would be to keep things as simple as possible - row-summation is the easier approach and induces less overhead in your code. You should remember that there is generally energy dissipation associated with using the lumped mass matrix regardless of the lumping method, so some refinement of timestep/element size might be needed depending on how well you need to conserve energy.

I've been using such a setup for solving dynamics problems with MPM (Material Point Method) and it works well Dirichlet Boundaries (i.e. a ball bouncing off of a rigid domain boundary).