Coarsening working in serial

src/GridapFixes.jl

@@ -45,4 +45,160 @@ function Base.map(::typeof(Gridap.Arrays.testitem),
   a1=Vector{eltype(eltype(a[1]))}(undef,size(a2,1))
   a1.=zero(Gridap.Arrays.testitem(a1))
   (a1,a2)
+end
+
+
+"""
+  Given a map from new cells to old cells, computes the inverse map.
+  In the general case (refinement+coarsening), the n2o map is a Table, 
+  but the algorithm is optimized for Vectors (refinement only).
+
+  IMPORTANT NOTE: This can go to Gridap.Adaptivity as far proper tests are written.
+
+"""
+function Gridap.Adaptivity.get_o2n_faces_map(ncell_to_ocell::Gridap.Arrays.Table{T}) where {T<:Integer}
+  println(ncell_to_ocell.data)
+  nC = maximum(ncell_to_ocell.data)
+  ptrs = fill(0,nC+1)
+  for ccell in ncell_to_ocell.data
+    ptrs[ccell+1] += 1
+  end
+  Gridap.Arrays.length_to_ptrs!(ptrs)
+  data = Vector{Int}(undef,ptrs[end]-1)
+  for fcell=1:length(ncell_to_ocell.ptrs)-1
+    for j=ncell_to_ocell.ptrs[fcell]:ncell_to_ocell.ptrs[fcell+1]-1
+      ccell=ncell_to_ocell.data[j]
+      data[ptrs[ccell]]=fcell
+      ptrs[ccell]+=1
+    end
+  end
+  Gridap.Arrays.rewind_ptrs!(ptrs)
+  ocell_to_ncell=Gridap.Arrays.Table(data,ptrs)
+  println(ocell_to_ncell.ptrs)
+  println(ocell_to_ncell.data)
+  ocell_to_ncell
+end
+
+"""
+  Given a AdaptivityGlue and a CellField defined on the parent(old) mesh, 
+  returns an equivalent CellField on the child(new) mesh.
+"""
+function Gridap.Adaptivity.change_domain_o2n(
+     f_old,
+     old_trian::Triangulation{Dc},
+     new_trian::AdaptedTriangulation,
+     glue::AdaptivityGlue) where {Dc}
+
+  oglue = get_glue(old_trian,Val(Dc))
+  nglue = get_glue(new_trian,Val(Dc))
+
+  Gridap.Helpers.@notimplementedif num_point_dims(old_trian) != Dc
+  Gridap.Helpers.@notimplementedif isa(nglue,Nothing)
+
+  if (num_cells(old_trian) != 0)
+    n2o_faces_map=glue.n2o_faces_map[end]
+    new_coarsened_cells_ids=findall(i->n2o_faces_map.ptrs[i+1]-n2o_faces_map.ptrs[i]>1,
+                                    1:length(n2o_faces_map.ptrs)-1)
+    # If mixed refinement/coarsening, then f_c2f is a Table
+    f_old_data=Gridap.CellData.get_data(f_old)
+    f_c2f=Gridap.Adaptivity.c2f_reindex(f_old_data,glue)
+    new_rrules = Gridap.Adaptivity.get_new_cell_refinement_rules(glue)
+    field_array=lazy_map(OldToNewField, f_c2f, new_rrules, glue.n2o_cell_to_child_id)
+    return Gridap.CellData.similar_cell_field(f_old,field_array,new_trian,ReferenceDomain())
+  else
+    f_new = Fill(Gridap.Fields.ConstantField(0.0),num_cells(new_trian))
+    return CellData.similar_cell_field(f_old,f_new,new_trian,ReferenceDomain())
+  end 
+end 
+
+function Gridap.Adaptivity.get_new_cell_refinement_rules(g::AdaptivityGlue{<:Gridap.Adaptivity.MixedGlue})
+  old_rrules = g.refinement_rules
+  n2o_faces_map = g.n2o_faces_map[end]
+  @assert isa(n2o_faces_map,Gridap.Arrays.Table)
+  lazy_map(Reindex(old_rrules), [n2o_faces_map.data[n2o_faces_map.ptrs[i]] for i=1:length(n2o_faces_map.ptrs)-1])
+end
+
+
+"""
+Given a domain and a non-overlapping refined cover, a `FineToCoarseField`
+is a `Field` defined in the domain and constructed by a set of fields defined on 
+the subparts of the covering partition.
+The refined cover is represented by a `RefinementRule`. 
+"""
+
+abstract type NewFieldType end;
+struct CoarsenedNewFieldType <: NewFieldType end;
+struct RefinedOrUntouchedNewFieldType <: NewFieldType end;   
+
+# Unfortunately, I cannot use traits in OldToNewField as its 
+# usage results in conversion errors when leveraging it with 
+# the lazy_map infrastructure
+struct OldToNewField <: Gridap.Fields.Field
+  new_field_type::NewFieldType
+  fine_to_coarse_field
+  refined_or_untouched_field
+  function OldToNewField(a::NewFieldType,b::Gridap.Fields.Field,c::Gridap.Fields.Field)
+    new(a,b,c)
+  end 
+end
+
+function OldToNewField(old_fields::AbstractArray{<:Gridap.Fields.Field},rrule::RefinementRule,child_id::Integer)
+  if length(old_fields)==1
+    cell_map = get_cell_map(rrule)[child_id]
+    old_field=old_fields[1]
+    fine_to_coarse_field=Gridap.Adaptivity.FineToCoarseField(
+          [old_field for i=1:Gridap.Adaptivity.num_subcells(rrule)],
+          rrule)
+    refined_or_untouched_field=old_field∘cell_map
+    OldToNewField(RefinedOrUntouchedNewFieldType(),fine_to_coarse_field,refined_or_untouched_field)
+  else 
+    Gridap.Helpers.@check length(old_fields) == Gridap.Adaptivity.num_subcells(rrule)
+    Gridap.Helpers.@check child_id==-1
+    cell_map = get_cell_map(rrule)[1]
+    fine_to_coarse_field=Gridap.Adaptivity.FineToCoarseField(old_fields,rrule)
+    refined_or_untouched_field=old_fields[1]∘cell_map 
+    OldToNewField(CoarsenedNewFieldType(),fine_to_coarse_field,refined_or_untouched_field)
+  end
+end
+
+function OldToNewField(old_field::Gridap.Fields.Field,rrule::RefinementRule,child_id::Integer)
+  Gridap.Helpers.@notimplemented
+end
+
+# # Necessary for distributed meshes, where not all children of a coarse cell may belong to the processor. 
+# function FineToCoarseField(fine_fields::AbstractArray{<:Field},rrule::RefinementRule,child_ids::AbstractArray{<:Integer})
+#   fields = Vector{Field}(undef,num_subcells(rrule))
+#   fields = fill!(fields,ConstantField(0.0))
+#   for (k,id) in enumerate(child_ids)
+#     fields[id] = fine_fields[k]
+#   end
+#   return FineToCoarseField(fields,rrule)
+# end
+
+function Gridap.Fields.return_cache(a::OldToNewField,x::AbstractArray{<:Point})
+  f2c_cache=Gridap.Fields.return_cache(a.fine_to_coarse_field,x)
+  rou_cache=Gridap.Fields.return_cache(a.refined_or_untouched_field,x)
+  return (f2c_cache,rou_cache)
+end
+
+function Gridap.Fields.evaluate!(cache,a::OldToNewField,x::AbstractArray{<:Point})
+  if isa(a.new_field_type,CoarsenedNewFieldType)
+    f2c_cache,rou_cache = cache
+    Gridap.Fields.evaluate!(f2c_cache,a.fine_to_coarse_field,x)
+  else
+    @assert isa(a.new_field_type,RefinedOrUntouchedNewFieldType)
+    f2c_cache,rou_cache = cache
+    Gridap.Fields.evaluate!(rou_cache,a.refined_or_untouched_field,x)
+  end  
+end
+
+# Fast evaluation of FineToCoarseFields: 
+# Points are pre-classified into the children cells, which allows for the search to be 
+# skipped entirely. 
+function Gridap.Fields.return_cache(a::OldToNewField,x::AbstractArray{<:Point},child_ids::AbstractArray{<:Integer})
+  Gridap.Helpers.@notimplemented
+end
+
+function Gridap.Fields.evaluate!(cache,a::OldToNewField,x::AbstractArray{<:Point},child_ids::AbstractArray{<:Integer})
+  Gridap.Helpers.@notimplemented
 end

src/OctreeDistributedDiscreteModels.jl

@@ -933,6 +933,7 @@ function _process_owned_cells_fine_to_coarse_model_glue(cmodel::DiscreteModel{Dc
   end
 
   function _setup_fine_to_coarse_faces_map_table(Dc,flags,num_o_c_cells,num_f_cells)
+    println("PPP: $(flags)")
     num_children = get_num_children(Val{Dc})
     fine_to_coarse_faces_map_ptrs = Vector{Int}(undef,num_f_cells+1)
     fine_to_coarse_faces_map_ptrs[1]=1
@@ -965,13 +966,13 @@ function _process_owned_cells_fine_to_coarse_model_glue(cmodel::DiscreteModel{Dc
     while cell <= num_o_c_cells
       if (flags[cell]==refine_flag)
         for child = 1:num_children
-          fine_to_coarse_faces_map_data[c+child-1] = cell
+          fine_to_coarse_faces_map_data[fine_to_coarse_faces_map_ptrs[c]+child-1] = cell
           fcell_to_child_id[c+child-1] = child
         end 
         c = c + num_children
         cell=cell+1
       elseif (flags[cell]==nothing_flag)
-        fine_to_coarse_faces_map_data[c]=cell
+        fine_to_coarse_faces_map_data[fine_to_coarse_faces_map_ptrs[c]]=cell
         fcell_to_child_id[c]=1
         c=c+1
         cell=cell+1
@@ -980,17 +981,22 @@ function _process_owned_cells_fine_to_coarse_model_glue(cmodel::DiscreteModel{Dc
         cell_fwd,coarsen=_move_fwd_and_check_if_all_children_coarsened(flags,num_o_c_cells,cell,num_children)
         if coarsen
           fcell_to_child_id[c]=-1
-          cell=cell+num_children
+          for child = 1:num_children
+            fine_to_coarse_faces_map_data[fine_to_coarse_faces_map_ptrs[c]+child-1] = cell
+            cell=cell+1
+          end 
+          c=c+1
         else 
           for j=cell:cell_fwd-1 
-            fine_to_coarse_faces_map_data[c]=j
+            fine_to_coarse_faces_map_data[fine_to_coarse_faces_map_ptrs[c]]=j
             fcell_to_child_id[c]=1
             c=c+1
             cell=cell+1
           end
         end
       end
     end
+    println("XXX: ", fine_to_coarse_faces_map_data)
     Gridap.Arrays.Table(fine_to_coarse_faces_map_data, fine_to_coarse_faces_map_ptrs),fcell_to_child_id
   end
 

test.jl

@@ -167,6 +167,57 @@ function test_coarsen(ranks,dmodel)
         flags
     end
     fmodel,glue=refine(dmodel,ref_coarse_flags);
+
+    order=1
+    reffe=ReferenceFE(lagrangian,Float64,order)
+    VH=FESpace(dmodel,reffe,conformity=:H1;dirichlet_tags="boundary")
+    UH=TrialFESpace(VH,u)
+
+    Vh=FESpace(fmodel,reffe,conformity=:H1;dirichlet_tags="boundary")
+    Uh=TrialFESpace(Vh,u)
+    ΩH  = Triangulation(dmodel)
+    dΩH = Measure(ΩH,degree)
+
+    aH(u,v) = ∫( ∇(v)⊙∇(u) )*dΩH
+    bH(v) = ∫(v*f)*dΩH
+
+    op = AffineFEOperator(aH,bH,UH,VH)
+    uH = solve(op)
+    e = u - uH
+
+    # # Compute errors
+    el2 = sqrt(sum( ∫( e*e )*dΩH ))
+    eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩH ))
+
+    tol=1e-8
+    @assert el2 < tol
+    @assert eh1 < tol
+
+
+    Ωh  = Triangulation(fmodel)
+    dΩh = Measure(Ωh,degree)
+
+    ah(u,v) = ∫( ∇(v)⊙∇(u) )*dΩh
+    bh(v) = ∫(v*f)*dΩh
+
+    op = AffineFEOperator(ah,bh,Uh,Vh)
+    uh = solve(op)
+    e = u - uh
+
+    # # Compute errors
+    el2 = sqrt(sum( ∫( e*e )*dΩh ))
+    eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩh ))
+
+    tol=1e-8
+    @assert el2 < tol
+    @assert eh1 < tol
+
+    # prolongation via interpolation
+    uHh=interpolate(uH,Uh)
+    e = uh - uHh
+    el2 = sqrt(sum( ∫( e*e )*dΩh ))
+    tol=1e-8
+    @assert el2 < tol
 end
 test_coarsen(ranks,dmodel);