Assembly Refactoring

[fverdugo]

Road map discussed with @oriolcg on 2021-02-1

1. Use AbstractUnitRage to describe dof ids in a FESpace since it will provide metadata needed e.g. for block arrays or partitioned arrays.

# Now
num_free_dofs(U::FESpace) = @abstractmethod
num_dirichlet_dofs(U::SingleFieldFESpace)  = @abstractmethod

# Update to
# Returns AbstractUnitRange
# In Single-field Base.OneTo
# In Multi-field MultiLevelBlockedUnitRange
get_free_dof_ids(U::FESpace) = @abstractmethod

# Returns Base.OneTo
get_dirichlet_dof_ids(U::SingleFieldFESpace) = @abstractmethod

num_free_dofs(U::FESpace) = length(get_free_dof_ids(U))
num_dirichlet_dofs(U::SingleFieldFESpace) = length(get_dirichlet_dof_ids(U))

2. Rename get_free_values -> get_free_dof_values and get_dirichlet_values -> get_dirichlet_dof_values to be more consistent with the new methods get_free_dof_ids and get_dirichlet_dof_ids.

3. Decouple FESpace and Assembler. This will remove current assumptions that preclude to fix issues like Multi-field in mixed-dimensional FE spaces #483

# Now,
get_test(a::Assembler) = @abstractmethod
get_trial(a::Assembler) = @abstractmethod

# Update to
#  -> AbstractUnitRange
get_rows(a::Assembler) = @abstractmethod
get_cols(a::Assembler) = @abstractmethod

num_rows(a::Assembler) = length(get_rows(a))
num_cols(a::Assembler) = length(get_cols(a))

Base.axes(a::Assembler) = (get_rows(a),get_cols(a))
Base.size(a::Assembler) = (num_rows(a),num_cols(a))

# Remove all usages of FE spaces in assemblers.

# Now,
function allocate_vector(a::SparseMatrixAssembler,vecdata)
  n = num_free_dofs(get_test(a))
  n = num_rows(a)
  allocate_vector(get_vector_type(a),n)
end

# Update to (Idem for all related cases )
function allocate_vector(a::SparseMatrixAssembler,vecdata)
  allocate_vector(get_vector_type(a),get_rows(a))
end

# Now
matdata = (cell_to_mat,cell_to_bgcell,cell_to_bgcell)

# Update to (Idem for similar cases)
# cell_to_mat already with constraints
# Idem for cell_to_row_ids and cell_to_col_ids
matdata = (cell_to_mat,cell_to_row_ids,cell_to_col_ids)


# Now
function collect_cell_matrix(mat_contributions)
  @abstractmethod
end

# Update to (Idem for similar cases)
function collect_cell_matrix(
  trial::FESpace,test::FESpace,mat_contributions)
  @abstractmethod
end

# Example:
function collect_cell_matrix(
  trial::FESpace,test::FESpace,a::DomainContribution)
  w = []
  r = []
  c = []
  for trian in get_domains(a)
    cell_mat = get_contribution(a,trian)
    @assert eltype(cell_mat) <: AbstractMatrix
    cell_to_bgcell = get_cell_to_bgcell(trian)
    cell_mat_c = attach_constraints_cols(trial,cell_mat,cell_to_bgcell)
    cell_mat_rc = attach_constraints_rows(test,cell_mat_c,cell_to_bgcell)
    rows = get_cell_dof_ids(test,cell_to_bgcell)
    cols = get_cell_dof_ids(trial,cell_to_bgcell)
    push!(w,cell_mat_rc)
    push!(r,rows)
    push!(c,cols)
  end
  (w,r,c)
end

4. Compress cell contributions. Needed e.g. to alleviate memory consumption in GridapEmbedded (cc @ericneiva @pmartorell ). Done in branch https://github.com/gridap/Gridap.jl/tree/compressed_assembly (We keep this uncoupled from the assembly refactoring of this issue for the moment)

# New method for in Triangulation
# Default implementation
function compress(cell_to_mat,trian::Triangulation)
  cell_to_mat, get_cell_to_bgcell(trian)
end

# Other triangulation types will specialize this to compress the `cell_to_mat` array when needed.  E.g.
function compress(cell_to_mat,trian::GridapEmbedded.SubCellTriangulation)
  # TODO return arrays indexed with a new cell id `ccell`, aka "compressed" cell
  ccell_to_mat, ccell_to_bgcell
end
 
# The new method will be called when collecting cell contributions eg:
function collect_cell_matrix(
  trial::FESpace,test::FESpace,a::DomainContribution)
  w = []
  r = []
  c = []
  for trian in get_domains(a)
    cell_mat = get_contribution(a,trian)
    @assert eltype(cell_mat) <: AbstractMatrix
    ccell_to_mat, ccell_to_bgcell = compress(cell_mat,trian)
    ccell_mat_c = attach_constraints_cols(trial,ccell_mat,ccell_to_bgcell)
    ccell_mat_rc = attach_constraints_rows(test,ccell_mat_c,ccell_to_bgcell)
    rows = get_cell_dof_ids(test,ccell_to_bgcell)
    cols = get_cell_dof_ids(trial,ccell_to_bgcell)
    push!(w,ccell_mat_rc)
    push!(r,rows)
    push!(c,cols)
  end
  (w,r,c)
end

• Re activate check about trial FE space (move it to AffineFEOperator / FEOperator)

[oriolcg]

I did step 3 in commit fcfd1c2. The AssemblersTests and SparseMatrixAssemblersTests are passing, but now we have an issue on the FEOperators. The problem is that the FEOperators get the trial and test FE spaces from the assembler via get_trial(assembler), which is not there anymore.

@fverdugo , what would be a good solution for that? Should we add the FESpaces as members of the FEOperators? Or just pass them as arguments where needed?

[fverdugo]

Should we add the FESpaces as members of the FEOperators?

Yes. It makes sense to have them as members.

[oriolcg]

Refactoring complete in 9322918. Should we do a PR and then merge with step 4 or merge before PR? @fverdugo @ericneiva @pmartorell

[fverdugo]

I would keep these two refactors separated for the moment. The work by @ericneiva @pmartorell can be readily included in a patch of v0.15. Your refactoring, would lead to 0.16.

[oriolcg]

Step 5 completed in c1f9a9b

[santiagobadia]

@oriolcg

Rename get_free_values -> get_free_dof_values and get_dirichlet_values -> get_dirichlet_dof_values to be more consistent with the new methods get_free_dof_ids and get_dirichlet_dof_ids.

I see that you have renamed the free values but not the fixed values. Now the notation is inconsistent. Can you fix it?

Thanks!

[fverdugo]

Some deprecation warnings would be useful here! So that people still using the old names know which new ones should use.

[fverdugo]

Closed via #606