Expanding adaptivity to non-BodyFittedTriangulation triangulations

NEWS.md

@@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Implemented `CompositeQuadrature`, a quadrature for a cell that has been refined using a `RefinementRule`. Since PR [#838](https://github.com/gridap/Gridap.jl/pull/838).
 - Implemented simple refinement strategies for Cartesian discrete models in 2&3D as well as Unstructured discrete models in 2D. The latter is implemented by red-green refinement. Since PR [#838](https://github.com/gridap/Gridap.jl/pull/838).
 - Added optimization when calling `unique` for a `CompressedArray`. Since PR [#838](https://github.com/gridap/Gridap.jl/pull/838).
+- Added support for changing domain between adapted triangulations in cases where the target triangulation is a view, a `BoundaryTriangulation` or a `SkeletonTriangulation`. Since PR [#868](https://github.com/gridap/Gridap.jl/pull/868).
 
 ### Fixed
 - Using broadcasting through in `ODESolver` vector operations. Since PR [#860](https://github.com/gridap/Gridap.jl/pull/860)

src/Adaptivity/AdaptedTriangulations.jl

@@ -112,14 +112,9 @@ for (stype,ttype) in [(:AdaptedTriangulation,:AdaptedTriangulation),(:AdaptedTri
         return is_change_possible($sstrian,$tttrian)
       end
 
-      if isa($sstrian,BodyFittedTriangulation) && isa($tttrian,BodyFittedTriangulation)
-        return is_related(strian,ttrian)
-      end
-
-      # Support for Boundary, Skeletton, ... triangulations is not yet implemented
-      # unless they are based on the same background model.
-      @notimplemented
-      return false
+      strian_is_cell_wise = (num_cell_dims(strian) == num_point_dims(strian))
+      trians_are_related  = is_related(strian,ttrian)
+      return strian_is_cell_wise && trians_are_related
     end
   end
   @eval begin
@@ -131,46 +126,49 @@ for (stype,ttype) in [(:AdaptedTriangulation,:AdaptedTriangulation),(:AdaptedTri
       if (get_background_model(strian) === get_background_model(ttrian))
         return best_target($sstrian,$tttrian)
       end
+
+      ttrian_is_cell_wise = (num_cell_dims(ttrian) == num_point_dims(ttrian))
+      strian_is_cell_wise = (num_cell_dims(strian) == num_point_dims(strian))
 
-      if isa($sstrian,BodyFittedTriangulation) && isa($tttrian,BodyFittedTriangulation)
+      if strian_is_cell_wise && ttrian_is_cell_wise
+        # Choose the child mesh, both ways are possible
         is_child(ttrian,strian) ? (return ttrian) : (return strian)
+      elseif strian_is_cell_wise
+        return strian
       end
-
-      # Support for Boundary, Skeletton, ... triangulations is not yet implemented
-      # unless they are based on the same background model.
+
       @notimplemented
       return nothing
     end
   end
 end
 
 function CellData.change_domain(a::CellField,strian::Triangulation,::ReferenceDomain,ttrian::AdaptedTriangulation,::ReferenceDomain)
-  if strian === ttrian
-    return a
-  end
   @check is_change_possible(strian,ttrian)
 
   if (get_background_model(strian) === get_background_model(ttrian))
-    return change_domain(a,strian,ReferenceDomain(),ttrian.trian,ReferenceDomain())
+    b = change_domain(a,strian,ReferenceDomain(),ttrian.trian,ReferenceDomain())
+    return CellData.similar_cell_field(b,get_data(b),ttrian,ReferenceDomain())
   end
 
-  return change_domain_o2n(a,ttrian)
+  return change_domain_o2n(a,strian,ttrian)
 end
 
 function CellData.change_domain(a::CellField,strian::AdaptedTriangulation,::ReferenceDomain,ttrian::Triangulation,::ReferenceDomain)
-  if strian === ttrian
-    return a
-  end
   @check is_change_possible(strian,ttrian)
 
   if (get_background_model(strian) === get_background_model(ttrian))
     return change_domain(a,strian.trian,ReferenceDomain(),ttrian,ReferenceDomain())
   end
 
-  return change_domain_n2o(a,ttrian)
+  return change_domain_n2o(a,strian,ttrian)
 end
 
 function CellData.change_domain(a::CellField,strian::AdaptedTriangulation,::ReferenceDomain,ttrian::AdaptedTriangulation,::ReferenceDomain)
+  if strian === ttrian
+    return a
+  end
+
   if is_child(strian,ttrian) # fine to coarse
     b = change_domain(a,strian,ReferenceDomain(),ttrian.trian,ReferenceDomain())
     return CellData.similar_cell_field(b,get_data(b),ttrian,ReferenceDomain())
@@ -191,12 +189,36 @@ for sdomain in [:ReferenceDomain,:PhysicalDomain]
   end
 end
 
-function Geometry.move_contributions(scell_to_val::AbstractArray, strian::AdaptedTriangulation, ttrian::Triangulation)
+function Geometry.move_contributions(scell_to_val::AbstractArray, strian::AdaptedTriangulation, ttrian::Triangulation{Dc}) where Dc
+  # Default Gridap 
+  if get_background_model(strian) === get_background_model(ttrian)
+    return move_contributions(scell_to_val,strian.trian,ttrian)
+  end
+
+  # Else 
   smodel = get_adapted_model(strian)
   @check get_parent(smodel) === get_background_model(ttrian)
 
-  tcell_to_val = move_contributions(scell_to_val,get_adaptivity_glue(smodel))
-  return tcell_to_val
+  return move_contributions_f2c(scell_to_val,strian,ttrian)
+end
+
+function move_contributions_f2c(fine_tface_to_val::AbstractArray, ftrian::AdaptedTriangulation, ctrian::Triangulation{Dc}) where Dc
+  glue  = get_adaptivity_glue(ftrian)
+  fglue = get_glue(ftrian,Val(Dc))
+  cglue = get_glue(ctrian,Val(Dc))
+
+  # Fine Triangulation -> Fine Model
+  fmodel = get_background_model(ftrian)
+  cell_ftrian = Triangulation(ReferenceFE{Dc},fmodel)
+  fine_mface_to_val = move_contributions(fine_tface_to_val,ftrian.trian,cell_ftrian)
+
+  # Fine model -> Coarse Model
+  coarse_mface_to_val = move_contributions(fine_mface_to_val,glue)
+
+  # Coarse Model -> Coarse Triangulation
+  coarse_tface_to_val = lazy_map(Reindex(coarse_mface_to_val),cglue.tface_to_mface)
+
+  return coarse_tface_to_val
 end
 
 function Geometry.move_contributions(scell_to_val::AbstractArray, glue::AdaptivityGlue)
@@ -206,72 +228,132 @@ function Geometry.move_contributions(scell_to_val::AbstractArray, glue::Adaptivi
   return tcell_to_val
 end
 
-function change_domain_o2n(f_old, new_trian::AdaptedTriangulation{Dc,Dp}) where {Dc,Dp}
-  @notimplementedif num_dims(get_triangulation(f_old)) != Dc
+# Change domain o2n/n2o
+
+function change_domain_o2n(f_old,old_trian::Triangulation,new_trian::AdaptedTriangulation)
   glue = get_adaptivity_glue(new_trian)
-  return change_domain_o2n(f_old,new_trian,glue)
+  return change_domain_o2n(f_old,old_trian,new_trian,glue)
 end
 
-function change_domain_n2o(f_new, old_trian::Triangulation{Dc,Dp}) where {Dc,Dp}
-  new_trian = get_triangulation(f_new)
-  @notimplementedif num_dims(new_trian) != Dc
-  @check isa(new_trian,AdaptedTriangulation)
+function change_domain_n2o(f_new,new_trian::AdaptedTriangulation,old_trian::Triangulation)
   glue = get_adaptivity_glue(new_trian)
-  return change_domain_n2o(f_new,old_trian,glue)
+  return change_domain_n2o(f_new,new_trian,old_trian,glue)
 end
 
 """
   Given a AdaptivityGlue and a CellField defined on the parent(old) mesh, 
   returns an equivalent CellField on the child(new) mesh.
 """
-function change_domain_o2n(f_old,new_trian::AdaptedTriangulation,glue::AdaptivityGlue)
+function change_domain_o2n(f_old,old_trian::Triangulation,new_trian::AdaptedTriangulation,glue::AdaptivityGlue)
   @notimplemented
 end
 
-function change_domain_o2n(f_coarse,ftrian::AdaptedTriangulation{Dc},glue::AdaptivityGlue{<:RefinementGlue,Dc}) where Dc
-  ctrian = get_triangulation(f_coarse)
-  @notimplementedif num_dims(ctrian) != Dc
+function change_domain_o2n(f_coarse,ctrian::Triangulation{Dc},ftrian::AdaptedTriangulation,glue::AdaptivityGlue{<:RefinementGlue}) where Dc
+  cglue = get_glue(ctrian,Val(Dc))
+  fglue = get_glue(ftrian,Val(Dc))
+
+  @notimplementedif num_point_dims(ctrian) != Dc
+  @notimplementedif isa(fglue,Nothing)
 
   if (num_cells(ctrian) != 0)
+    ### Old Triangulation -> Old Model
+    coarse_tface_to_field = CellData.get_data(f_coarse)
+    coarse_mface_to_field = extend(coarse_tface_to_field,cglue.mface_to_tface)
+
+    ### Old Model -> New Model
     # Coarse field but with fine indexing, i.e 
     #   f_f2c[i_fine] = f_coarse[coarse_parent(i_fine)]
-    f_f2c = c2f_reindex(f_coarse,glue)
+    f_f2c = c2f_reindex(coarse_mface_to_field,glue)
 
     # Fine to coarse coordinate map: x_coarse = Φ^(-1)(x_fine)
     ref_coord_map = get_n2o_reference_coordinate_map(glue)
 
-    # Final map: f_fine(x_fine) = f_f2c ∘ Φ^(-1)(x_fine) = f_coarse(x_coarse)
-    f_fine = lazy_map(∘,f_f2c,ref_coord_map)
-    return GenericCellField(f_fine,ftrian,ReferenceDomain())
+    # f_fine(x_fine) = f_f2c ∘ Φ^(-1)(x_fine) = f_coarse(x_coarse)
+    fine_mface_to_field = lazy_map(∘,f_f2c,ref_coord_map)
+
+    ### New Model -> New Triangulation
+    fine_tface_to_field = lazy_map(Reindex(fine_mface_to_field),fglue.tface_to_mface)
+    f_fine = lazy_map(Broadcasting(∘),fine_tface_to_field,fglue.tface_to_mface_map)
+
+    return CellData.similar_cell_field(f_coarse,f_fine,ftrian,ReferenceDomain())
   else
     f_fine = Fill(Gridap.Fields.ConstantField(0.0),num_cells(ftrian))
-    return GenericCellField(f_fine,ftrian,ReferenceDomain())
+    return CellData.similar_cell_field(f_coarse,f_fine,ftrian,ReferenceDomain())
   end
 end
 
 """
   Given a AdaptivityGlue and a CellField defined on the child(new) mesh, 
   returns an equivalent CellField on the parent(old) mesh.
 """
-function change_domain_n2o(f_new,old_trian::Triangulation,glue::AdaptivityGlue)
+function change_domain_n2o(f_new,new_trian::AdaptedTriangulation,old_trian::Triangulation,glue::AdaptivityGlue)
   @notimplemented
 end
 
-function change_domain_n2o(f_fine,ctrian::Triangulation{Dc},glue::AdaptivityGlue{<:RefinementGlue,Dc}) where Dc
-  @notimplementedif num_dims(ctrian) != Dc
-  msg = "Evaluating a fine CellField in the coarse mesh is costly! If you are using this feature 
-         to integrate, consider using a CompositeMeasure instead (see test/AdaptivityTests/GridTransferTests.jl)."
-  @warn msg
+function change_domain_n2o(f_fine,ftrian::AdaptedTriangulation{Dc},ctrian::Triangulation,glue::AdaptivityGlue{<:RefinementGlue}) where Dc
+  fglue = get_glue(ftrian,Val(Dc))
+  cglue = get_glue(ctrian,Val(Dc))
+
+  @notimplementedif num_point_dims(ftrian) != Dc
+  @notimplementedif isa(cglue,Nothing)
 
   if (num_cells(ctrian) != 0)
+    ### New Triangulation -> New Model
+    fine_tface_to_field = CellData.get_data(f_fine)
+    fine_mface_to_field = extend(fine_tface_to_field,fglue.mface_to_tface)
+
+    ### New Model -> Old Model
     # f_c2f[i_coarse] = [f_fine[i_fine_1], ..., f_fine[i_fine_nChildren]]
-    f_c2f = f2c_reindex(f_fine,glue)
+    f_c2f = f2c_reindex(fine_mface_to_field,glue)
 
     rrules   = get_old_cell_refinement_rules(glue)
-    f_coarse = lazy_map(FineToCoarseField,f_c2f,rrules)
-    return GenericCellField(f_coarse,ctrian,ReferenceDomain())
+    coarse_mface_to_field = lazy_map(FineToCoarseField,f_c2f,rrules)
+
+    ### 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))
+    return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain())
+  end
+end
+
+
+# Specialisation for Skeleton Pairs
+
+function change_domain_o2n(f_old,old_trian::Triangulation,new_trian::AdaptedTriangulation{Dc,Dp,<:SkeletonTriangulation}) where {Dc,Dp}
+  @check isa(CellData.get_data(f_old),Fill)
+  new_trian_plus = AdaptedTriangulation(new_trian.trian.plus,new_trian.adapted_model)
+  return change_domain_o2n(f_old,old_trian,new_trian_plus)
+end
+
+function change_domain_n2o(f_new,new_trian::AdaptedTriangulation,old_trian::SkeletonTriangulation)
+  @check isa(CellData.get_data(f_new),Fill)
+  return change_domain_n2o(f_new,new_trian,old_trian.plus)
+end
+
+function change_domain_o2n(f_old::CellData.CellFieldAt,old_trian::Triangulation,new_trian::AdaptedTriangulation{Dc,Dp,<:SkeletonTriangulation}) where {Dc,Dp}
+  if isa(f_old,CellData.CellFieldAt{:plus})
+    _new_trian = AdaptedTriangulation(new_trian.trian.plus,new_trian.adapted_model)
+  elseif isa(f_old,CellData.CellFieldAt{:minus})
+    _new_trian = AdaptedTriangulation(new_trian.trian.minus,new_trian.adapted_model)
   else
-    f_coarse = Fill(Gridap.Fields.ConstantField(0.0),num_cells(ftrian))
-    return GenericCellField(f_coarse,ctrian,ReferenceDomain())
+    @unreachable
   end
+  return change_domain_o2n(f_old,old_trian,_new_trian)
 end
+
+function change_domain_n2o(f_new::CellData.CellFieldAt,new_trian::AdaptedTriangulation,old_trian::SkeletonTriangulation)
+  if isa(f_new,CellData.CellFieldAt{:plus})
+    _old_trian = old_trian.plus
+  elseif isa(f_new,CellData.CellFieldAt{:minus})
+    _old_trian = old_trian.minus
+  else
+    @unreachable
+  end
+  return change_domain_n2o(f_new,new_trian,_old_trian)
+end
+
+

src/Adaptivity/Adaptivity.jl

@@ -47,7 +47,6 @@ include("AdaptivityGlues.jl")
 include("AdaptedDiscreteModels.jl")
 include("AdaptedTriangulations.jl")
 include("CompositeQuadratures.jl")
-
 include("EdgeBasedRefinement.jl")
 
 end # module