Getting global indices for element mass matrix per cell

[ruboerner]

TLDR; A grid with P cells is considered. How to obtain a representation of the element mass matrices arranged as a vector of sparse matrices, where each of these P matrices would essentially be the result of a matrix assembly based on a grid with a parameterization that is defined by the i-th unit vector, i=1,...,P.

Problem statement

I am working on a geophysical parameter estimation problem. The parameter to be recovered is the spatial distribution of the electrical conductivity, which appears as a linear coefficient in Maxwell's curl-curl equation. After FE discretization, the action of the conductivity on the solution is materialized in the (N x N) mass matrix $\mathbf{M}$.

With $\mathbf{m}$ as P-vector of piece-wise constant cell conductivities, we get the following discretized PDE

$$\left(\mathbf{K} + i\omega \mathbf{M}(\mathbf{m})\right) \mathbf{u}(\mathbf{m}) = \mathbf{f}.$$

Sensitivity

To calculate the sensitivity of the physical field with respect to some (small) change in the conductivity, we need the partial derivative of the numerical solution $\mathbf{u}'(\mathbf{m})$ with respect to the P individual cell conductivities. Since the conductivity appears linearly in the weak form, we could form the element mass matrices (equivalent to an assignment of a piece-wise uniform parametrization with the value of 1) arranged as a very sparse matrix (or any other convenient form) for each of the P individual cells. This forms a sparse (N x N x P) 3-way tensor, or a P-vector of sparse N x N-arrays.

A further building block for the sensitivity calculation is the multiplication of the FE solution along the second dimension of this tensor (each "slice" of which is an N x N matrix) followed by arranging these P products into an N x P array resulting in

$$\mathbf{M}'(\mathbf{m}) \times_2 \mathbf{u}$$

(Sidenote: The multiplication with the P cells along the third dimension and subsequent summation over the scaled slices would be algebraically equivalent to the assembly of the mass matrix.)

It would be nice to have a construct that stores individual element mass matrices with entries spread out at the appropriate rows and columns just as they would finally arrive in a complete assembly.

(Sidenote 2: The tensor-representation and the associated tensor-times-vector routines have already been implemented and used in an earlier project, where an in-house FE library has been developed and implemented in Julia.)

The problem is: Although I can access the individual element matrices as well as the mapping between local and global DOF indices using get_cell_dof_ids(), I get lost with the negative entries for the dirichlet dofs. Further, there seems to be a different indexing scheme for free and dirichlet dofs.

A discretization with second order Lagrange elements and 18 triangular cells results in

l2g_map = get_cell_dof_ids(V, Ω) # local-to-global mapping

julia> l2g_map
18-element Table{Int32, Vector{Int32}, Vector{Int32}}:
 [-1, -2, -5, -13, -14, 5]
 [-2, -5, 1, 5, 6, 7]
 [-2, -3, 1, -15, 6, 8]
 [-3, 1, 2, 8, 9, 10]
 [-3, -4, 2, -16, 9, 11]
 [-4, 2, -6, 11, -17, 12]
 [-5, 1, -7, 7, -18, 13]
 [1, -7, 3, 13, 14, 15]
 [1, 2, 3, 10, 14, 16]
 [2, 3, 4, 16, 17, 18]
 [2, -6, 4, 12, 17, 19]
 [-6, 4, -8, 19, -19, 20]
 [-7, 3, -9, 15, -20, 21]
 [3, -9, -10, 21, 22, -21]
 [3, 4, -10, 18, 22, 23]
 [4, -10, -11, 23, 24, -22]
 [4, -8, -11, 20, 24, 25]
 [-8, -11, -12, 25, -23, -24]

Where can I find - or how to implement - the appropriate mapping from this representation to one which is based on the global numbering of DOFs?

I would be extremely grateful for help.

Best regards,
Ralph

[JordiManyer]

It would be nice to have a construct that stores individual element mass matrices with entries spread out at the appropriate rows and columns just as they would finally arrive in a complete assembly.

It's basically the get_cell_dof_ids. The dirichlet dof ids (negative indices) do not get assembled in the matrix, they get assembled in the rhs vector. You can have a look at src/CellData/AttachDirichlet.jl. When assembling the matrix, the negative indices get ignored. So the cell dof ids really give you the global rows/cols of the matrix. Just not everything gets assembled.

[ruboerner]

Thank you very much for the quick response! It makes totally sense.

In case anybody might find it useful, here is an MWE with my (very quick and possibly naive) solution:

using Gridap
using Gridap.Geometry
using Gridap.FESpaces
using LinearAlgebra
using SparseArrays

"""
    assemble_mass_slice(
        local_to_global_ids::Vector{Int32},
        M_ref::Matrix{Float64},
        sz::Int) -> SparseMatrixCSC

Assemble an individual reference element matrix `M_ref` with 
respective entries according to their global indexing defined in 
`local_to_global_ids`. The function returns an (sz x sz) slice of a 
sparse (sz x sz x P) three-way tensor.
"""
function assemble_mass_slice(local_to_global_ids, M_ref, sz)
  nD = length(local_to_global_ids)
  M_ref = M_ref[:]
  Cij = zeros(Int, nD * nD, 2)
  Cij .= reduce(
      vcat,
      [
          [x[1] x[2]] for
          x in collect(Iterators.product(local_to_global_ids, local_to_global_ids))
      ],
  )
  mask = Bool.(collect(prod(Cij .> 0, dims = 2)))
  ind = findall(mask .== true)
  row_ind = [i[1] for i in ind]
  A = sparse(Cij[row_ind, 1], Cij[row_ind, 2], M_ref[row_ind], sz, sz)

  return A
end

model = UnstructuredDiscreteModel(
    CartesianDiscreteModel((0, 1, 0, 1), (3, 3))) |> simplexify
Ω = Triangulation(model)

order = 2
reffe = ReferenceFE(lagrangian, Float64, order)
V = TestFESpace(model,
    reffe;
    conformity=:H1,
    dirichlet_tags="boundary")

g(x) = -1.0
U = TrialFESpace(V, g)

num_dofs = num_free_dofs(V)

degree = 2 * order
dΩ = Measure(Ω, degree)

a(u, v) = ∫( v ⋅ u )*dΩ
M_ref = assemble_matrix(a, U, V)

du = get_trial_fe_basis(U)
dv = get_fe_basis(V)

data = collect_cell_matrix(U, V, a(du, dv))

elem_matrices = data[1][1]
indices = data[2][1]
n_cells = length(indices)

A = Vector{SparseMatrixCSC{Float64, Int64}}(undef, n_cells)

for i in 1:n_cells
    A[i] = assemble_mass_slice(indices[i], elem_matrices[i], num_dofs)
end

M_manual = sum(A, dims=1)[1]

@show norm(M_manual .- M_ref)

[JordiManyer]

No worries, glad it helped!