I am trying to use NGSolve to realize something similar to Section 4.2 of the paper by Alexander Linke and Christian Merdon, where a mass matrix of the BDM-interpolated discontinuous FE space is added to the bilinear operator of Navier Stokes equations. May I ask if there is any way to directly add to the BilinearForm formula a mass matrix with trial and test bases interpolated/projected into another space, e.g. BDM? Right now I can separately get the mass matrix with bases interpolated by embedding into the target space and back into the original discontinuous FE space, but I do not know how to add this separate mass matrix to the bilinear form to get the operator inverse. Could you provide some insight into this? Thank you for your time!