Merge pull request #886 from gridap/refinement-rule-subface-glue

Set of improvements/fixes in Gridap.Adaptivity module

Author: amartinhuertas

NEWS.md

@@ -10,11 +10,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 - Jacobi polynomial bases. Since PR [#896](https://github.com/gridap/Gridap.jl/pull/896).
-
 - Replaced newest vertex bisection mesh adaptation in
   `src/Geometry/NewestVertexBisection.jl` with appropriate changes to
-  `src/Adaptivity/EdgeBasedRefinement.jl`. see PR
+  `src/Adaptivity/EdgeBasedRefinement.jl`. Since PR
   [#901](https://github.com/gridap/Gridap.jl/pull/901).
+- When refining `DiscreteModels`, the `FaceLabeling` of the resulting `AdaptedDiscreteModel` will now correctly inhering the tags of the parent model. This has been made possible by the addition of the method `get_d_to_face_to_parent_face`. Since PR[#886](https://github.com/gridap/Gridap.jl/pull/886).
+- Added support for mixed adaptivity (i.e coarsening and refining), as well as non-conforming adaptivity. Since PR[#886](https://github.com/gridap/Gridap.jl/pull/886).
 
 ### Changed
 
@@ -26,9 +27,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Fixed the method `get_normal_vector` for `AdaptedTriangulation`. The method `get_facet_normal`
   was using default, it's now using the spetialized implementation for the underlying triangulation type.
   Since PR [#884](https://github.com/gridap/Gridap.jl/pull/884).
-
 - Fixed `cell_dof_ids` for the case of vectorial `ConstantFESpace`. Since PR [#888](https://github.com/gridap/Gridap.jl/pull/888)
 - Fixed generation of Modal C0 bases for Julia 1.9. Since PR [#918](https://github.com/gridap/Gridap.jl/pull/918).
+- Fixed some edge cases for `change_domain` between `AdaptedTriangulations` where inneficient coordinate transformations would be applied between physical and reference domains. Since PR[#886](https://github.com/gridap/Gridap.jl/pull/886).
+- Fixed: Domain limits can now be of any type (notably, floats) when refining `CartesianDiscreteModels`. Since PR[#886](https://github.com/gridap/Gridap.jl/pull/886).
 
 ## [0.17.17] - 2023-02-24
 

src/Adaptivity/AdaptedDiscreteModels.jl

@@ -103,24 +103,31 @@ function refine(model::CartesianDiscreteModel{Dc}, cell_partition::Tuple) where
   desc = Geometry.get_cartesian_descriptor(model)
   nC   = desc.partition
 
-  # Refined model
-  domain    = _get_cartesian_domain(desc)
-  model_ref = CartesianDiscreteModel(domain,cell_partition.*nC)
-
-  # Glue
+  # Refinement Glue
   f2c_cell_map, fcell_to_child_id = _create_cartesian_f2c_maps(nC,cell_partition)
-  faces_map      = [(d==Dc) ? f2c_cell_map : Int[] for d in 0:Dc]
-  reffe          = LagrangianRefFE(Float64,first(get_polytopes(model)),1)
-  rrules         = RefinementRule(reffe,cell_partition)
+  faces_map = [(d==Dc) ? f2c_cell_map : Int[] for d in 0:Dc]
+  reffe     = LagrangianRefFE(Float64,first(get_polytopes(model)),1)
+  rrules    = RefinementRule(reffe,cell_partition)
   glue = AdaptivityGlue(faces_map,fcell_to_child_id,rrules)
 
+  # Refined model
+  domain     = _get_cartesian_domain(desc)
+  _model_ref = CartesianDiscreteModel(domain,cell_partition.*nC)
+
+  # Propagate face labels
+  coarse_labels = get_face_labeling(model)
+  coarse_topo   = get_grid_topology(model)
+  fine_topo     = get_grid_topology(_model_ref)
+  fine_labels   = _refine_face_labeling(coarse_labels,glue,coarse_topo,fine_topo)
+
+  model_ref = CartesianDiscreteModel(get_grid(_model_ref),fine_topo,fine_labels)
   return AdaptedDiscreteModel(model_ref,model,glue)
 end
 
 function _get_cartesian_domain(desc::CartesianDescriptor{D}) where D
   origin = desc.origin
   corner = origin + VectorValue(desc.sizes .* desc.partition)
-  domain = Vector{Int}(undef,2*D)
+  domain = Vector{eltype(origin)}(undef,2*D)
   for d in 1:D
     domain[d*2-1] = origin[d]
     domain[d*2]   = corner[d]
@@ -159,3 +166,10 @@ end
     return f2c_map, child_map
   end)
 end
+
+function get_d_to_fface_to_cface(model::AdaptedDiscreteModel)
+  ftopo = get_grid_topology(get_model(model))
+  ctopo = get_grid_topology(get_parent(model))
+  glue  = get_adaptivity_glue(model)
+  return get_d_to_fface_to_cface(glue,ctopo,ftopo)
+end

src/Adaptivity/AdaptedTriangulations.jl

@@ -183,8 +183,14 @@ end
 
 for sdomain in [:ReferenceDomain,:PhysicalDomain]
   for (stype,ttype) in [(:AdaptedTriangulation,:AdaptedTriangulation),(:AdaptedTriangulation,:Triangulation),(:Triangulation,:AdaptedTriangulation)]
+    sstrian = (stype==:AdaptedTriangulation) ? :(strian.trian) : :(strian)
+    tttrian = (ttype==:AdaptedTriangulation) ? :(ttrian.trian) : :(ttrian)
     @eval begin
-      function CellData.change_domain(a::CellField,strian::$stype,::$sdomain,ttrian::$ttype,::PhysicalDomain)
+      function CellData.change_domain(a::CellField,strian::$stype,sd::$sdomain,ttrian::$ttype,::PhysicalDomain)
+        if (get_background_model(strian) === get_background_model(ttrian))
+          return change_domain(a,$sstrian,sd,$tttrian,PhysicalDomain())
+        end
+
         a_ref  = change_domain(a,ReferenceDomain())
         atrian = change_domain(a_ref,strian,ReferenceDomain(),ttrian,ReferenceDomain())
         return change_domain(atrian,PhysicalDomain())
@@ -281,11 +287,36 @@ function change_domain_o2n(f_coarse,ctrian::Triangulation{Dc},ftrian::AdaptedTri
 
     return CellData.similar_cell_field(f_coarse,f_fine,ftrian,ReferenceDomain())
   else
-    f_fine = Fill(Gridap.Fields.ConstantField(0.0),num_cells(ftrian))
+    f_fine = Fill(Fields.ConstantField(0.0),num_cells(ftrian))
     return CellData.similar_cell_field(f_coarse,f_fine,ftrian,ReferenceDomain())
   end
 end
 
+function change_domain_o2n(
+  f_old,
+  old_trian::Triangulation{Dc},
+  new_trian::AdaptedTriangulation,
+  glue::AdaptivityGlue{<:MixedGlue}) where {Dc}
+
+  oglue = get_glue(old_trian,Val(Dc))
+  nglue = get_glue(new_trian,Val(Dc))
+
+  @notimplementedif num_point_dims(old_trian) != Dc
+  @notimplementedif isa(nglue,Nothing)
+
+  if (num_cells(old_trian) != 0)
+    # If mixed refinement/coarsening, then f_c2f is a Table
+    f_old_data  = CellData.get_data(f_old)
+    f_c2f       = c2f_reindex(f_old_data,glue)
+    new_rrules  = get_new_cell_refinement_rules(glue)
+    field_array = lazy_map(OldToNewField, f_c2f, new_rrules, glue.n2o_cell_to_child_id)
+    return CellData.similar_cell_field(f_old,field_array,new_trian,ReferenceDomain())
+  else
+    f_new = Fill(Fields.ConstantField(0.0),num_cells(new_trian))
+    return CellData.similar_cell_field(f_old,f_new,new_trian,ReferenceDomain())
+  end 
+end
+
 """
   Given a AdaptivityGlue and a CellField defined on the child(new) mesh, 
   returns an equivalent CellField on the parent(old) mesh.
@@ -310,16 +341,17 @@ function change_domain_n2o(f_fine,ftrian::AdaptedTriangulation{Dc},ctrian::Trian
     # f_c2f[i_coarse] = [f_fine[i_fine_1], ..., f_fine[i_fine_nChildren]]
     f_c2f = f2c_reindex(fine_mface_to_field,glue)
 
-    rrules   = get_old_cell_refinement_rules(glue)
-    coarse_mface_to_field = lazy_map(FineToCoarseField,f_c2f,rrules)
+    child_ids = f2c_reindex(glue.n2o_cell_to_child_id,glue)
+    rrules    = get_old_cell_refinement_rules(glue)
+    coarse_mface_to_field = lazy_map(FineToCoarseField,f_c2f,rrules,child_ids)
 
     ### Old Model -> Old Triangulation
     coarse_tface_to_field = lazy_map(Reindex(coarse_mface_to_field),cglue.tface_to_mface)
     f_coarse = lazy_map(Broadcasting(∘),coarse_tface_to_field,cglue.tface_to_mface_map)
 
     return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain())
   else
-    f_coarse = Fill(Gridap.Fields.ConstantField(0.0),num_cells(fcoarse))
+    f_coarse = Fill(Fields.ConstantField(0.0),num_cells(fcoarse))
     return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain())
   end
 end

src/Adaptivity/Adaptivity.jl

@@ -42,6 +42,7 @@ export change_domain, move_contributions
 
 include("RefinementRules.jl")
 include("FineToCoarseFields.jl")
+include("OldToNewFields.jl")
 include("FineToCoarseReferenceFEs.jl")
 include("AdaptivityGlues.jl")
 include("AdaptedDiscreteModels.jl")

src/Adaptivity/AdaptivityGlues.jl

@@ -22,13 +22,20 @@ struct AdaptivityGlue{GT,Dc,A,B,C,D,E} <: GridapType
   is_refined           :: E
 
   function AdaptivityGlue(n2o_faces_map::Vector{<:Union{AbstractVector{<:Integer},Table{<:Integer}}},
-                          n2o_cell_to_child_id::AbstractVector{<:Integer},
+                          n2o_cell_to_child_id::Union{AbstractVector{<:Integer},Table{<:Integer}},
                           refinement_rules::AbstractVector{<:RefinementRule})
     Dc = length(n2o_faces_map)-1
     is_refined    = select_refined_cells(n2o_faces_map[Dc+1])
     o2n_faces_map = get_o2n_faces_map(n2o_faces_map[Dc+1])
 
     GT = all(is_refined) ? RefinementGlue : MixedGlue
+    if isa(GT(),RefinementGlue)
+      @assert isa(n2o_faces_map,Vector{<:AbstractVector{<:Integer}})
+      @assert isa(n2o_cell_to_child_id,AbstractVector{<:Integer})
+    else
+      @assert isa(n2o_faces_map,Vector{<:Table{<:Integer}})
+      @assert isa(n2o_cell_to_child_id,Table{<:Integer})
+    end
 
     A = typeof(n2o_faces_map)
     B = typeof(n2o_cell_to_child_id)
@@ -67,7 +74,26 @@ end
   but the algorithm is optimized for Vectors (refinement only).
 """
 function get_o2n_faces_map(ncell_to_ocell::Table{T}) where {T<:Integer}
-  @notimplemented
+  nC = maximum(ncell_to_ocell.data)
+  
+  ptrs = fill(0,nC+1)
+  for ccell in ncell_to_ocell.data
+    ptrs[ccell+1] += 1
+  end
+  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
+  Arrays.rewind_ptrs!(ptrs)
+
+  ocell_to_ncell = Table(data,ptrs)
+  return ocell_to_ncell
 end
 
 function get_o2n_faces_map(ncell_to_ocell::Vector{T}) where {T<:Integer}
@@ -95,12 +121,19 @@ function get_o2n_faces_map(ncell_to_ocell::Vector{T}) where {T<:Integer}
   return ocell_to_ncell
 end
 
-function get_new_cell_refinement_rules(g::AdaptivityGlue)
+function get_new_cell_refinement_rules(g::AdaptivityGlue{<:RefinementGlue})
   old_rrules = g.refinement_rules
   n2o_faces_map = g.n2o_faces_map[end]
   return lazy_map(Reindex(old_rrules),n2o_faces_map)
 end
 
+function get_new_cell_refinement_rules(g::AdaptivityGlue{<:MixedGlue})
+  old_rrules = g.refinement_rules
+  n2o_faces_map = g.n2o_faces_map[end]
+  new_idx = lazy_map(Reindex(n2o_faces_map.data),n2o_faces_map.ptrs[1:end-1])
+  return lazy_map(Reindex(old_rrules), new_idx)
+end
+
 function get_old_cell_refinement_rules(g::AdaptivityGlue)
   return g.refinement_rules
 end
@@ -129,3 +162,104 @@ function _reindex(data,idx::Vector)
   m = Reindex(data)
   return lazy_map(m,idx)
 end
+
+
+# New to old face glues 
+
+function get_d_to_fface_to_cface(::AdaptivityGlue,::GridTopology,::GridTopology)
+  @notimplemented
+end
+
+"
+ For each child/fine face, returns the parent/coarse face containing it. The parent 
+ face might have higher dimension. 
+
+ Returns two arrays: 
+  - [dimension][fine face gid] -> coarse parent face gid
+  - [dimension][fine face gid] -> coarse parent face dimension
+"
+function get_d_to_fface_to_cface(glue::AdaptivityGlue{<:RefinementGlue},
+                                 ctopo::GridTopology{Dc},
+                                 ftopo::GridTopology{Dc}) where Dc
+
+  # Local data for each coarse cell, at the RefinementRule level
+  rrules = Adaptivity.get_old_cell_refinement_rules(glue)
+  ccell_to_d_to_faces = lazy_map(rr->map(d->Geometry.get_faces(get_grid_topology(rr.ref_grid),Dc,d),0:Dc),rrules)
+  ccell_to_d_to_fface_to_parent_face = lazy_map(get_d_to_face_to_parent_face,rrules)
+
+  # Global data, concerning the complete meshes
+  ccell_to_fcell = glue.o2n_faces_map
+  d_to_ccell_to_cface = map(d->Geometry.get_faces(ctopo,Dc,d),0:Dc)
+  d_to_fcell_to_fface = map(d->Geometry.get_faces(ftopo,Dc,d),0:Dc)
+
+  d_to_fface_to_cface = [fill(Int32(0),num_faces(ftopo,d)) for d in 0:Dc]
+  d_to_fface_to_cface_dim = [fill(Int32(0),num_faces(ftopo,d)) for d in 0:Dc]
+
+  # For each coarse cell
+  for ccell in 1:num_cells(ctopo)
+    local_d_to_fface_to_parent_face,
+      local_d_to_fface_to_parent_dim = ccell_to_d_to_fface_to_parent_face[ccell]
+    # For each fine subcell:
+    # child_id -> Local Id of the fine cell within the refinement rule (ccell)
+    for (child_id,fcell) in enumerate(ccell_to_fcell[ccell])
+      # For each fine face on the fine subcell: 
+      # d     -> Dimension of the fine face
+      # iF    -> Local Id of the fine face within the fine cell
+      # fface -> Global Id of the fine face 
+      for d in 0:Dc
+        for (iF,fface) in enumerate(d_to_fcell_to_fface[d+1][fcell])
+          # Local Id of the fine face within the refinement rule
+          fface_child_id = ccell_to_d_to_faces[ccell][d+1][child_id][iF]
+          # Local Id of the coarse parent face within the coarse cell
+          parent    = local_d_to_fface_to_parent_face[d+1][fface_child_id]
+
+          # Global Id of the coarse parent face, and it's dimension
+          cface_dim = local_d_to_fface_to_parent_dim[d+1][fface_child_id]
+          cface     = d_to_ccell_to_cface[cface_dim+1][ccell][parent]
+          d_to_fface_to_cface[d+1][fface] = cface
+          d_to_fface_to_cface_dim[d+1][fface] = cface_dim
+        end
+      end
+    end
+  end
+
+  return (d_to_fface_to_cface,d_to_fface_to_cface_dim)
+end
+
+# FaceLabeling refinement
+
+function _refine_face_labeling(coarse_labeling::FaceLabeling,
+                               glue  :: AdaptivityGlue,
+                               ctopo :: GridTopology,
+                               ftopo :: GridTopology)
+  d_to_fface_to_cface,
+    d_to_fface_to_cface_dim = get_d_to_fface_to_cface(glue,ctopo,ftopo)
+
+  return _refine_face_labeling(coarse_labeling,d_to_fface_to_cface,d_to_fface_to_cface_dim)
+end
+
+function _refine_face_labeling(coarse_labeling::FaceLabeling,
+                               d_to_fface_to_cface,
+                               d_to_fface_to_cface_dim)
+  tag_to_name = copy(coarse_labeling.tag_to_name)
+  tag_to_entities = copy(coarse_labeling.tag_to_entities)
+  
+  Dc = num_dims(coarse_labeling)
+  d_to_dface_to_entity = Vector{Vector{Int32}}(undef,Dc+1)
+  for d in 0:Dc
+    nF = length(d_to_fface_to_cface[d+1])
+    dface_to_entity = Vector{Int32}(undef,nF)
+  
+    for fface in 1:nF
+      cface = d_to_fface_to_cface[d+1][fface]
+      cface_dim = d_to_fface_to_cface_dim[d+1][fface]
+
+      cface_entity = coarse_labeling.d_to_dface_to_entity[cface_dim+1][cface]
+      dface_to_entity[fface] = cface_entity
+    end
+  
+    d_to_dface_to_entity[d+1] = dface_to_entity
+  end
+  
+  return Geometry.FaceLabeling(d_to_dface_to_entity,tag_to_entities,tag_to_name)  
+end