Bulk penalty stab div

[amartinhuertas]

@amboschman ... Mostly creating this PR so that I can add comments to the code

not to be merged YET .... WIP

[amartinhuertas]

@amboschman ... I have taken a careful look at the lines below this one till the call to revised_setup_aggregate_to_cut_cells. I think I also got what the issue/bug with your original setup_aggregate_to_cut_cells was. At the end of the day, you need to generate an array with the root_cells and an array with the cut_cells in the very same order you have them listed in aggregate_to_cells

In any case, I think there is a better (more general and more efficient) way to do the very same thing. I copy and paste the code that implements this alternative way below:

# Set up bg_cell_type
# -1: exterior 
#  0: cut 
#  1: interior, non root 
# >1: interior root 

# TO-DO: pack this for loop in a function `set_up_cell_type`
cell_type = zeros(Int64,length(aggregates))
for (pos,val) in enumerate(aggregates)
   if val==0
      cell_type[pos]=-1
   else 
      cell_type[val]+=1
   end
end

cell_is_exterior(cell)=cell_type[cell]==-1
cell_is_cut(cell)=cell_type[cell]==0
cell_is_interior_root(cell)=cell_type[cell]>1
cell_is_interior_non_root(cell)=cell_type[cell]==1
 
function flatten(array_of_arrays)
   vcat(array_of_arrays...)
end

function restrict_to_cell_type(aggregate_to_cells, cell_type_filter)
    ret=Vector{Int}[]
    for agg in aggregate_to_cells
      new_agg=Int64[]
      for cell in agg
         if cell_type_filter(cell)
            push!(new_agg, cell)
         end 
      end 
      push!(ret, new_agg)
    end
    ret
end 

With this machinery, you can compute root_cells, cut_cells and aggregate_to_cut_cells in a unified way as:

root_cells = flatten(restrict_to_cell_type(aggregate_to_cells,cell_is_interior_root))
aggregate_to_cut_cells = restrict_to_cell_type(aggregate_to_cells,cell_is_cut))
cut_cells   = flatten(aggregate_to_cut_cells)

Let me know if that makes sense.