"Augment Lagrangian fails to solve & Wrong lumped mass matrix" for higher-order elements

[Huangzizhou]

Describe the bug
The BCLagrangianForm uses the lumped mass matrix as the weights for Augmented Lagrangian. For linear elements, the lumped mass matrix is diagonally dominant, so the weights are positive; for higher-order elements, the mass matrix is not necessarily diagonally dominant, though still PSD, so the weights can be zero or even negative and cannot be used.

The same issue happens to the transient simulations if the lumped mass is used. We need a proper way to lump the mass matrix for higher-order elements.

[zfergus]

Copying this discussion from Slack for reference:

Lumping for higher-order elements cannot be done by summing rows (this works only for P1): the resulting lumped matrix may be not PD, even if the full mass is PD.

A simple solution is diagonal scaling [1] which is discarding all off-diagonal elements in each element mass matrix, and then rescaling the diagonal ones so that they sum up to the element area/volume, then assembling the complete matrix from per-element in the standard way.

An even simpler solution (which looses more accuracy, but may be enough for AL, but not for time stepping) is to assign to each node the sum of the areas/volumes of all incident elements, scaled by 1/(number of nodes per element), e.g. 6 for quadratic triangle, and 10 for a quadratic tet.

[1] Hinton, E., T. Rock, and O. C. Zienkiewicz. "A note on mass lumping and related processes in the finite element method." Earthquake Engineering & Structural Dynamics 4.3 (1976): 245-249