Function for assembling L2MM and invL2MM in one shot

[amartinhuertas]

Hi @davelee2804 @cbladwell,

find below a code snippet to assemble the (possibly block diagonal) L2 mass matrix and its inverse "in one shot". I think this is an efficient way to go. Note that:

1. The code avoids extracting the diagonal blocks from an already assembled global sparse matrix. This is very inefficient. (scalar sparse matrix storage formats are not well prepared for such kind of operations).
2. The code collects the contribution of each cell to the global matrix into a regular vector (i.e., no longer lazy) using collect, i.e., explicitly storing all cell matrices at once (i.e., the matrix diagonal blocks). Usually we cannot afford this. However, in this particular case, the amount of memory required is equivalent to the one storing all global matrix entries. As a result, we avoid computing the cell matrices twice (once for L2MM and another one for invL2MM).
3. We use similar in order to get a sparse matrix with the same sparsity pattern as L2MM but with undefined entries yet. This is far more efficient than assembling a matrix from scratch. (We have already commented this in other contexts).

Let me know if everything makes sense ....

using Gridap
using GridapGeosciences

order=1
degree=4

model = CubedSphereDiscreteModel(1; radius=rₑ)

Ω = Triangulation(model)
dΩ = Measure(Ω, degree)

reffe_lgn = ReferenceFE(lagrangian, Float64, order)
Q = FESpace(model, reffe_lgn; conformity=:L2)
P = TrialFESpace(Q)

function assemble_L2MM_invL2MM(U,V,dΩ)
  a(a,b) = ∫(a⋅b)dΩ
  v = get_fe_basis(V)
  u = get_trial_fe_basis(U)
  assemblytuple    = Gridap.FESpaces.collect_cell_matrix(U,V,a(u,v))
  cell_matrix_MM   = collect(assemblytuple[1][1]) # This result is no longer a LazyArray
  newassemblytuple = ([cell_matrix_MM],assemblytuple[2],assemblytuple[3])
  a=SparseMatrixAssembler(U,V)
  L2MM=assemble_matrix(a,newassemblytuple)
  for i in eachindex(cell_matrix_MM)
     cell_matrix_MM[i]=inv(cell_matrix_MM[i])
  end
  invL2MM=similar(L2MM)
  Gridap.FESpaces.assemble_matrix!(invL2MM, a, newassemblytuple)
  L2MM, invL2MM
end

A,B=assemble_L2MM_invL2MM(P,Q,dΩ)
println(A*B)

[davelee2804]

Thanks @amartinhuertas , I've added this to the new shallow_water_imex branch

[amartinhuertas]

Thanks @amartinhuertas , I've added this to the new shallow_water_imex branch

Ok, cool, closing the issue, and opening PR ...