* canvas code now outputs the aggregates to vtk. Bounding boxes for the

  aggregrates seem to be correctly computer.

* Fixed assembly of lhs of the local L2 projection into the bounding box space.
  Now we integrate on the domain defined as the union of the aggregated
  cells, not on the domain of the whole bounding box as before.

bulk_ghost_penalty_canvas.jl

@@ -1,5 +1,6 @@
 using Gridap
 using GridapEmbedded
+using FillArrays
 
 # Manufactured solution
 order = 1
@@ -12,7 +13,6 @@ ud(x) = u(x)
 R = 0.2
 geom = disk(R, x0=Point(0.5,0.5))
 
-
 # Setup background model
 n=10
 partition = (n,n)
@@ -111,27 +111,61 @@ function setup_aggregates_bounding_box_model(bgmodel, aggregate_to_cells)
 end 
 
 aggregate_to_cells=setup_aggregate_to_cells(aggregates)
-aaggregates_bounding_box_model=
+aggregates_bounding_box_model=
        setup_aggregates_bounding_box_model(bgmodel,aggregate_to_cells)
 
+colors = color_aggregates(aggregates,bgmodel)
+writevtk(Triangulation(bgmodel),"trian",celldata=["cellin"=>aggregates,"color"=>colors])       
 writevtk(aggregates_bounding_box_model, "bb_model")
 writevtk(bgmodel, "bg_model")
 
-# Compute LHS of L2 projection 
-order=1
-degree=2*order
-# To use Q instead of P! 
-reffe=ReferenceFE(lagrangian,Float64,order)
-# We need a discontinuous space to represent the L2 projection
-Vbb=FESpace(aggregates_bounding_box_model,reffe,conformity=:L2)
-Ubb=TrialFESpace(Vbb)
-Ωbb=Triangulation(aggregates_bounding_box_model)
-dΩbb=Measure(Ωbb,degree)
-ubb=get_trial_fe_basis(Ubb)
-vbb=get_fe_basis(Vbb)
-lhs=get_array(∫(vbb*ubb)dΩbb)
+"""
+  Changes the domain of a trial/test basis defined on 
+  the reference space of bounding boxes to the reference 
+  space of the agg cells
+
+  TO-DO: in the future, for system of PDEs (MultiField) we should 
+         also take care of blocks (BlockMap)
+"""
+function change_domain_bb_to_agg_cells(basis_bb, 
+                                       ref_agg_cell_to_ref_bb_map, 
+                                       Ωagg_cells,
+                                       agg_cells_to_aggregate)
+    @assert num_cells(Ωagg_cells)==length(ref_agg_cell_to_ref_bb_map)
+    @assert Gridap.CellData.DomainStyle(basis_bb)==ReferenceDomain()
+    bb_basis_style = Gridap.FESpaces.BasisStyle(basis_bb)  
+    bb_basis_array = Gridap.CellData.get_data(basis_bb)
+    if (bb_basis_style==Gridap.FESpaces.TrialBasis())
+        # Remove transpose map; we will add it later
+        @assert isa(bb_basis_array,Gridap.Arrays.LazyArray)
+        @assert isa(bb_basis_array.maps,Fill)
+        @assert isa(bb_basis_array.maps.value,typeof(transpose))
+        bb_basis_array=bb_basis_array.args[1]
+    end
+
+    bb_basis_array_to_Ωagg_cells_array = lazy_map(Reindex(bb_basis_array),agg_cells_to_aggregate)
+    bb_basis_array_to_Ωagg_cells_array = lazy_map(Broadcasting(∘),
+                                                  bb_basis_array_to_Ωagg_cells_array,
+                                                  ref_agg_cell_to_ref_bb_map)
+    if (bb_basis_style==Gridap.FESpaces.TrialBasis())
+        # Add transpose 
+        bb_basis_array_to_Ωagg_cells_array=lazy_map(transpose, bb_basis_array_to_Ωagg_cells_array)
+    end    
+
+    Gridap.CellData.GenericCellField(bb_basis_array_to_Ωagg_cells_array,
+                                     Ωagg_cells,
+                                     ReferenceDomain())
+end
+
+
+# Set up objects required to compute both LHS and RHS of the L2 projection
 
-# Compute RHS of L2 projection
+# Set up global spaces 
+Ωhact = Triangulation(cutdisk,ACTIVE)
+order = 1
+reffe = ReferenceFE(lagrangian,Float64,order)
+Vstd  = TestFESpace(Ωhact,reffe,conformity=:L2)
+Ustd  = TrialFESpace(Vstd)
 
 # Generate an array with the global IDs of the cells 
 # that belong to an aggregrate. From now on, we will 
@@ -150,6 +184,7 @@ ref_agg_cell_to_agg_cell_map=get_cell_map(Ωagg_cells)
 
 # Generate an array that given the local ID of an "agg_cell"
 # returns the ID of the aggregate to which it belongs
+# (i.e., flattened version of aggregrate_to_cells)
 agg_cells_to_aggregate=Vector{Int}()
 for (i,cells) in enumerate(aggregate_to_cells)
     for _ in cells 
@@ -163,23 +198,70 @@ bb_to_ref_bb_agg_cells=lazy_map(Reindex(bb_to_ref_bb),agg_cells_to_aggregate)
 ref_agg_cell_to_ref_bb_map=
     lazy_map(Broadcasting(∘),bb_to_ref_bb_agg_cells,ref_agg_cell_to_agg_cell_map)
 
-# Change domain of vbb from Ωbb to Ωagg_cells
-vbb_array=Gridap.CellData.get_data(vbb)
-vbb_array_to_Ωagg_cells_array = lazy_map(Reindex(vbb_array),agg_cells_to_aggregate)
-vbb_array_to_Ωagg_cells_array = lazy_map(Broadcasting(∘),
-                                         vbb_array_to_Ωagg_cells_array,
-                                         ref_agg_cell_to_ref_bb_map)
-vbb_Ωagg_cells=Gridap.CellData.GenericCellField(vbb_array_to_Ωagg_cells_array,
-                                                Ωagg_cells,
-                                                ReferenceDomain())                                              
+
+# Compute LHS of L2 projection 
+degree=2*order
+dΩagg_cells = Measure(Ωagg_cells,degree)
+order=1
+reffe=ReferenceFE(lagrangian,Float64,order) # Here we MUST use a Q space (not a P space!)
+Vbb=FESpace(aggregates_bounding_box_model,reffe,conformity=:L2) # We need a DG space to represent the L2 projection
+Ubb=TrialFESpace(Vbb)
+ubb=get_trial_fe_basis(Ubb)
+vbb=get_fe_basis(Vbb)
+
+# Change domain of vbb (test) from Ωbb to Ωagg_cells
+vbb_Ωagg_cells=change_domain_bb_to_agg_cells(vbb,
+                                             ref_agg_cell_to_ref_bb_map,
+                                             Ωagg_cells,
+                                             agg_cells_to_aggregate)
+# Change domain of ubb (trial) from Ωbb to Ωagg_cells                                             
+ubb_Ωagg_cells=change_domain_bb_to_agg_cells(ubb,
+                                             ref_agg_cell_to_ref_bb_map,
+                                             Ωagg_cells,
+                                             agg_cells_to_aggregate)
+
+# Compute contributions to LHS of L2 projection
+agg_cells_to_lhs_contribs=get_array(∫(vbb_Ωagg_cells*ubb_Ωagg_cells)dΩagg_cells)
+
+aggregate_to_local_cells=copy(aggregate_to_cells)
+current_local_cell=1
+for (i,cells) in enumerate(aggregate_to_local_cells)
+    for j in 1:length(cells)
+        cells[j]=current_local_cell
+        current_local_cell+=1
+    end
+end
+
+# TO-DO: Better name?
+struct AssembleLhsMap{A} <: Gridap.Fields.Map
+    agg_cells_lhs_contribs::A
+end 
+
+function Gridap.Fields.return_cache(m::AssembleLhsMap,cells)
+    cache_unassembled_lhs=array_cache(m.agg_cells_lhs_contribs)
+    evaluate_result=Gridap.Arrays.CachedArray(eltype(eltype(m.agg_cells_lhs_contribs)),2)
+    cache_unassembled_lhs,evaluate_result
+end
+
+function Gridap.Fields.evaluate!(cache,m::AssembleLhsMap,cells)
+    cache_unassembled_lhs,result=cache
+    contrib = getindex!(cache_unassembled_lhs,m.agg_cells_lhs_contribs,1)
+    Gridap.Arrays.setsize!(result,(size(contrib,1),size(contrib,2)))
+    result.array .= 0.0
+    for (i,cell) in enumerate(cells)
+        contrib = getindex!(cache_unassembled_lhs,m.agg_cells_lhs_contribs,cell)
+        result.array .= contrib
+    end
+    result.array
+end
+
+# Finally assemble LHS contributions
+ass_lhs_map=AssembleLhsMap(agg_cells_to_lhs_contribs)
+lhs=lazy_map(ass_lhs_map,aggregate_to_local_cells)
 
 # Compute contributions to the RHS of the L2 projection
-Ωhact = Triangulation(cutdisk,ACTIVE)
-Vstd  = TestFESpace(Ωhact,reffe,conformity=:L2)
-Ustd  = TrialFESpace(Vstd)
 du    = get_trial_fe_basis(Ustd)
 Ωagg_cell_dof_ids = get_cell_dof_ids(Ustd,Ωagg_cells)
-dΩagg_cells = Measure(Ωagg_cells,degree)
 agg_cells_rhs_contribs=get_array(∫(vbb_Ωagg_cells*du)dΩagg_cells)
 
 # TO-DO: Better name?
@@ -230,15 +312,6 @@ function Gridap.Fields.evaluate!(cache,m::AssembleRhsMap,cells)
   result.array
 end
 
-aggregate_to_local_cells=copy(aggregate_to_cells)
-current_local_cell=1
-for (i,cells) in enumerate(aggregate_to_local_cells)
-    for j in 1:length(cells)
-        cells[j]=current_local_cell
-        current_local_cell+=1
-    end
-end
-
 ass_rhs_map=AssembleRhsMap(Ωagg_cell_dof_ids,agg_cells_rhs_contribs)
 rhs=lazy_map(ass_rhs_map,aggregate_to_local_cells)