Hi all,
I am doing Reduced-Order Modeling (ROM) with legacy FEniCS for a nonlinear PDE. I would like to assemble specific entries (or sets of entries) of a matrix rather than the full matrix, with the aim of extracting the set of assembled matrix entries for subsequent calculations (for use in a DEIM-type ROM approach, for those who know). I have a working approach that does this in a roundabout way by modifying the integration measure to integrate over only those cells which support the matrix entries that I’m interested in.
While that approach works and gives me correct results, it also results in a lot of “bycatch”, since all other basis functions that are supported on the selected cells are evaluated as well. I suspect that those redundant evaluations are slowing down my calculations significantly (relative to the regular matrix assembly approach). I’m aware of the RBniCS package, which seemingly has an efficient implementation of DEIM where (I assume) individual (sets of) matrix entries are being assembled, though I haven’t been able to figure out from RBniCS’ source code how this is achieved in an efficient way. Does FEniCS even allow for such an approach without extensive modifications to the assembly routines?
The way I do it in RBniCS is by creating a submesh, see
If you want to do it yourself, you may be better off using recent contributions by @cdaversin, since that has the capability of building submeshes (called MeshView) and assembling on them.
Thank you for mentioning this file! Am I correct in saying that you start with a MeshFunction and create a submesh (and associated FunctionSpace) that contains only the cells you are interested in? Have you found any advantages of this approach over directly using the MeshFunction to modify the integration measure of a form? To me (as a relative novice) it seems like both approaches should result in similar operations (where we assemble over a subset of all cells in the mesh), but I’m not familiar with performance advantages/drawbacks of either approach.
Using a submesh/meshview and then assemble gives a smaller linear system.
If you use the function space on the full mesh, your matrix/vector will contain Zero entries wherever you have not integrated over, and thus your system might not be invertible.
That’s a good point. In my use case having that larger matrix is desirable. I’m using tIGAr to do Isogeometric Analysis (IGA) with (legacy) FEniCS. Having the full matrix available allows any zero entries from the finite element system to be propagated to the IGA system. The IGA system does not have to be invertible, since in this reduced-order modeling application we don’t solve it directly but instead process it in some way.