Code structure: Assembly

System matrix assembly should be done with the minimum requirements of information possible. Specific properties of each object used during assembling should be the responsibility of the object itself, not the assembly algorithm.

Assembly algorithm for bilinear forms

• mixed_space: MixedSpace containing n FESpaces.
• S_local: n x n local matrices corresponding to the terms relating two pairs of spaces in mixed_space, for example, for the mixed Poisson equation n = 2 (a 2-form $q^{(2)}$ and a 3-form $\phi^{(3)}$). We will have a local dense matrix blocks $S^{11} = [\langle\epsilon^{(2)}_i, \epsilon^{(2)}_j\rangle]$, $S^{12} = [-\langle\nabla\cdot\epsilon^{(2)}_i, \epsilon^{(3)}_j\rangle]$, $S^{21} = [\langle\epsilon^{(3)}_i, \nabla\cdot\epsilon^{(2)}_j\rangle]$, $S^{22} = [0]$.
• dof_offset: each FESpace in a MixedFEMSpace has a dof_offset specifying where in the system matrix its degrees of freedom should start. This vector contains the offsets for each of the spaces in mixed_space.
• n_spaces: the number of spaces making up mixed_space.

for element_idx in element_list
    # Compute all cross inner products between the active bases
    # (\varphi_{i}, \varphi_{j})_{\Omega_{element_idx}}

    # Use the element assembler to compute the local evaluation of the bilinear form
    active_bases, n_spaces, dof_offset, S_local = element_assembly(element_idx, mixed_space, data_fields)

    # Populate the index vector and value vector of the sparse matrix representing the bilinear form
    for space_l in 1:n_spaces, space_r in 1:n_spaces
        if isnothing(S_local(space_l, space_r))
            continue
        end
        for (base_l_idx, base_l) in enumerate(active_bases), (base_r_idx, base_r) in enumerate(active_bases)
            S_i_idx[data_idx] = base_l + space_l * dof_offset[space_l]  # for a mixed space we need to offset the indexing of the dofs
            S_j_idx[data_idx] = base_r + space_r * dof_offset[space_r]
            S_values[data_idx] = S_local[base_l_idx, base_r_idx, space_l, space_r]  # need to check best indexing, but we need to index the different local matrices that come up from the interaction between the different spaces
        end
    end
end

# Assemble the sparse matrix
S = sparse(S_i_idx, S_j_idx, S_values)

Open topics

• This algorithm does not explicitly include the mapping. How do we wish to treat that?

  • Do we implement that in the Assembly function?
  • Or do we construct wedge and star operators (as proposed in the algorithm) and use that? If we add these two operators then we can use the algorithm. Is this efficient?

• We need different, but similar, algorithms for different assemblies

  • Mass matrix assembly
  • Stiffness matrix assembly
  • Forcing term assembly
  • Nonlinear terms assembly (and other different terms)

• The basis needs to include its rank (k-form) so that we can compute wedges between different k-forms.