Merge pull request #932 from gridap/adaptivity-docs

Adaptivity docs

docs/make.jl

@@ -19,6 +19,7 @@ pages = [
   "Gridap.Visualization" => "Visualization.md",
   "Gridap.FESpaces" => "FESpaces.md",
   "Gridap.MultiField" => "MultiField.md",
+  "Gridap.Adaptivity" => "Adaptivity.md",
   "Developper notes" => Any[
     "dev-notes/block-assemblers.md",
   ],

docs/src/Adaptivity.md

@@ -0,0 +1,112 @@
+# Gridap.Adaptivity
+
+```@meta
+CurrentModule = Gridap.Adaptivity
+```
+
+The adaptivity module provides a framework to work with adapted (refined/coarsened/mixed) meshes.
+
+It provides
+
+- A generic interface to represent adapted meshes and a set of tools to work with Finite Element spaces defined on them. In particular, moving `CellFields` between parent and child meshes.
+- Particular implementations for conformally refining/coarsening 2D/3D meshes using several well-known strategies. In particular, Red-Green refinement and longest-edge bisection.
+
+## Interface
+
+The following types are defined in the module:
+
+```@docs
+RefinementRule
+AdaptivityGlue
+AdaptedDiscreteModel
+AdaptedTriangulation
+```
+
+The high-level interface is provided by the following methods:
+
+```@docs
+refine
+coarsen
+adapt
+```
+
+## Edge-Based refinement
+
+The module provides a `refine` method for `UnstructuredDiscreteModel`. The method takes a string `refinement_method`
+that determines the refinement strategy to be used. The following strategies are available:
+
+- `"red_green"` :: Red-Green refinement, default.
+- `"nvb"` :: Longest-edge bisection (only for meshes of TRIangles)
+
+Additionally, the method takes a kwarg `cells_to_refine` that determines which cells will be refined.
+Possible input types are:
+
+- `Nothing` :: All cells get refined.
+- `AbstractArray{<:Bool}` of size `num_cells(model)` :: Only cells such that `cells_to_refine[iC] == true` get refined.
+- `AbstractArray{<:Integer}` :: Cells for which `gid ∈ cells_to_refine` get refined
+
+The algorithms try to respect the `cells_to_refine` input as much as possible, but some additional cells
+might get refined in order to guarantee that the mesh remains conforming.
+
+```julia
+  function refine(model::UnstructuredDiscreteModel;refinement_method="red_green",kwargs...)
+    [...]
+  end
+```
+
+## CartesianDiscreteModel refining
+
+The module provides a `refine` method for `CartesianDiscreteModel`. The method takes a `Tuple` of size `Dc`
+(the dimension of the model cells) that will determine how many times cells will be refined in
+each direction. For example, for a 2D model, `refine(model,(2,3))` will refine each QUAD cell into
+a 2x3 grid of cells.
+
+```julia
+  function refine(model::CartesianDiscreteModel{Dc}, cell_partition::Tuple) where Dc
+    [...]
+  end
+```
+
+## Notes for users
+
+Most of the tools provided by this module are showcased in the tests of the module itself, as well as the following tutorial (coming soon).
+
+However, we want to stress a couple of key performance-critical points:
+
+- The refining/coarsening routines are not optimized for performance. In particular, they are not parallelized.
+  If you require an optimized/parallel implementation, please consider leveraging specialised meshing libraries. For instance, we provide an implementation of `refine/coarsen` using P4est in the [GridapP4est.jl](https://github.com/gridap/GridapP4est.jl) library.
+
+- Although the toolbox allows you to evaluate `CellFields` defined on both fine/coarse meshes on their parent/children mesh, both directions of evaluation are not equivalent. As a user, you should always try to evaluate/integrate on the finest mesh for maximal performance. Evaluating a fine `CellField` on a coarse mesh relies on local tree searches, and is therefore a very expensive operation that should be avoided whenever possible.
+
+## Notes for developers
+
+### RefinementRule API
+
+Given a `RefinementRule`, the library provides a set of methods to compute the mappings between parent (coarse) face ids and child (fine) face ids (and vice-versa). The ids are local to the `RefinementRule`.
+
+```@docs
+get_d_to_face_to_child_faces
+get_d_to_face_to_parent_face
+get_face_subface_ldof_to_cell_ldof
+```
+
+### AdaptivityGlue API
+
+```@docs
+get_n2o_reference_coordinate_map
+get_old_cell_refinement_rules
+get_new_cell_refinement_rules
+get_d_to_fface_to_cface
+n2o_reindex
+o2n_reindex
+```
+
+### New-to-old field evaluations
+
+When a cell is refined, we need to be able to evaluate the fields defined on the children cells on the parent cell. To do so, we bundle the fields defined on the children cells into a new type of `Field` called `FineToCoarseField`. When evaluated on a `Point`, a `FineToCoarseField` will select the child cell that contains the `Point` and evaluate the mapped point on the corresponding child field.
+
+```@docs
+FineToCoarseField
+FineToCoarseDofBasis
+FineToCoarseRefFE
+```

docs/src/index.md

@@ -52,6 +52,6 @@ Pages = [
   "Visualization.md",
   "FESpaces.md",
   "MultiField.md",
+  "Adaptivity.md",
   ]
 ```
-

src/Adaptivity/AdaptedDiscreteModels.jl

@@ -1,8 +1,8 @@
 
 """
-  AdaptedDiscreteModel
 
-  `DiscreteModel` created by refining/coarsening another `DiscreteModel`. 
+  `DiscreteModel` created by refining/coarsening another `DiscreteModel`.
+  
   The refinement/coarsening hierarchy can be traced backwards by following the 
   `parent` pointer chain. This allows the transfer of dofs 
   between `FESpaces` defined on this model and its ancestors.
@@ -51,6 +51,11 @@ is_related(m1::DiscreteModel,m2::DiscreteModel) = is_child(m1,m2) || is_child(m2
 
 # Model Adaptation
 
+"""
+  function refine(model::DiscreteModel,args...;kwargs...) :: AdaptedDiscreteModel
+
+  Returns an `AdaptedDiscreteModel` that is the result of refining the given `DiscreteModel`.
+"""
 function refine(model::DiscreteModel,args...;kwargs...) :: AdaptedDiscreteModel
   @abstractmethod
 end
@@ -60,10 +65,21 @@ function refine(model::AdaptedDiscreteModel,args...;kwargs...)
   return AdaptedDiscreteModel(ref_model.model,model,ref_model.glue)
 end
 
+"""
+  function coarsen(model::DiscreteModel,args...;kwargs...) :: AdaptedDiscreteModel
+
+  Returns an `AdaptedDiscreteModel` that is the result of coarsening the given `DiscreteModel`.
+"""
 function coarsen(model::DiscreteModel,args...;kwargs...) :: AdaptedDiscreteModel
   @abstractmethod
 end
 
+"""
+  function adapt(model::DiscreteModel,args...;kwargs...) :: AdaptedDiscreteModel
+
+  Returns an `AdaptedDiscreteModel` that is the result of adapting (mixed coarsening and refining) 
+  the given `DiscreteModel`.
+"""
 function adapt(model::DiscreteModel,args...;kwargs...) :: AdaptedDiscreteModel
   @abstractmethod
 end

src/Adaptivity/AdaptedTriangulations.jl

@@ -1,13 +1,12 @@
 
 """
-  AdaptedTriangulation
 
   Triangulation produced from an AdaptedDiscreteModel.
   
   Contains: 
 
-  - adapted_model ~> AdaptedDiscreteModel for the triangulation.
-  - trian ~> Triangulation extracted from the background model, i.e get_model(adapted_model).
+  - adapted_model :: `AdaptedDiscreteModel` for the triangulation.
+  - trian :: `Triangulation` extracted from the background model, i.e `get_model(adapted_model)`.
 """
 struct AdaptedTriangulation{Dc,Dp,A<:Triangulation{Dc,Dp},B<:AdaptedDiscreteModel} <: Triangulation{Dc,Dp}
   trian::A
@@ -273,7 +272,7 @@ function change_domain_o2n(f_coarse,ctrian::Triangulation{Dc},ftrian::AdaptedTri
     ### 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(coarse_mface_to_field,glue)
+    f_f2c = o2n_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)
@@ -307,7 +306,7 @@ function change_domain_o2n(
   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)
+    f_c2f       = o2n_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())
@@ -339,9 +338,9 @@ function change_domain_n2o(f_fine,ftrian::AdaptedTriangulation{Dc},ctrian::Trian
 
     ### New Model -> Old Model
     # f_c2f[i_coarse] = [f_fine[i_fine_1], ..., f_fine[i_fine_nChildren]]
-    f_c2f = f2c_reindex(fine_mface_to_field,glue)
+    f_c2f = n2o_reindex(fine_mface_to_field,glue)
 
-    child_ids = f2c_reindex(glue.n2o_cell_to_child_id,glue)
+    child_ids = n2o_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)
 

src/Adaptivity/AdaptivityGlues.jl

@@ -4,15 +4,24 @@ struct RefinementGlue <: AdaptivityGlueType end
 struct MixedGlue <: AdaptivityGlueType end
 
 """
-Adaptivity glue between two nested triangulations:
-
-- `n2o_faces_map`       : Given a new face gid, returns 
-        A) if fine, the gid of the old face containing it.
-        B) if coarse, the gids of its children (in child order)
-- `n2o_cell_to_child_id`: Given a new cell gid, returns 
-        A) if fine, the local child id within the (old) coarse cell containing it.
-        B) if coarse, a list of local child ids of the (old) cells containing it.
-- `refinement_rules`    : RefinementRule used for each coarse cell.
+Glue containing the map between two nested triangulations. The contained datastructures will 
+depend on the type of glue. There are two types of `AdaptivityGlue`: 
+
+- `RefinementGlue` :: All cells in the new mesh are children of cells in the old mesh. I.e given 
+  a new cell, it is possible to find a single old cell containing it (the new cell might be exactly
+  the old cell if no refinement).
+- `MixedGlue` :: Some cells in the new mesh are children of cells in the old mesh, while others are
+  parents of cells in the old mesh. 
+
+Contains: 
+
+- `n2o_faces_map` :: Given a new face gid, returns
+  - if fine, the gid of the old face containing it.
+  - if coarse, the gids of its children (in child order)
+- `n2o_cell_to_child_id` :: Given a new cell gid, returns 
+  - if fine, the local child id within the (old) coarse cell containing it.
+  - if coarse, a list of local child ids of the (old) cells containing it.
+- `refinement_rules` :: Array conatining the `RefinementRule` used for each coarse cell.
 """
 struct AdaptivityGlue{GT,Dc,A,B,C,D,E} <: GridapType
   n2o_faces_map        :: A
@@ -62,7 +71,7 @@ select_refined_cells(n2o_cell_map::Vector) = Fill(true,length(n2o_cell_map))
 select_refined_cells(n2o_cell_map::Table) = map(x -> length(x) == 1, n2o_cell_map)
 
 """
-For each fine cell, returns Φ st. x_coarse = ϕ(x_fine)
+For each fine cell, returns the map Φ st. x_coarse = ϕ(x_fine)
 """
 function get_n2o_reference_coordinate_map(g::AdaptivityGlue{RefinementGlue})
   rrules    = get_new_cell_refinement_rules(g)
@@ -127,6 +136,10 @@ function get_o2n_faces_map(ncell_to_ocell::Vector{T}) where {T<:Integer}
   return ocell_to_ncell
 end
 
+"""
+  Given a `RefinementGlue`, returns an array containing the refinement rules for each
+  cell in the new mesh. 
+"""
 function get_new_cell_refinement_rules(g::AdaptivityGlue{<:RefinementGlue})
   old_rrules = g.refinement_rules
   n2o_faces_map = g.n2o_faces_map[end]
@@ -140,24 +153,38 @@ function get_new_cell_refinement_rules(g::AdaptivityGlue{<:MixedGlue})
   return lazy_map(Reindex(old_rrules), new_idx)
 end
 
+"""
+  Given a `RefinementGlue`, returns an array containing the refinement rules for each
+  cell in the old mesh. 
+"""
 function get_old_cell_refinement_rules(g::AdaptivityGlue)
   return g.refinement_rules
 end
 
 # Data re-indexing
 
-function f2c_reindex(fine_data,g::AdaptivityGlue)
-  ccell_to_fcell = g.o2n_faces_map
-  return _reindex(fine_data,ccell_to_fcell)
+"""
+  `function n2o_reindex(new_data,g::AdaptivityGlue) -> old_data`
+
+  Reindexes a cell-wise array from the new mesh to the old mesh.
+"""
+function n2o_reindex(new_data,g::AdaptivityGlue)
+  ocell_to_ncell = g.o2n_faces_map
+  return _reindex(new_data,ocell_to_ncell)
 end
 
-function c2f_reindex(coarse_data,g::AdaptivityGlue{GT,Dc}) where {GT,Dc}
-  fcell_to_ccell = g.n2o_faces_map[Dc+1]
-  return _reindex(coarse_data,fcell_to_ccell)
+"""
+  `function o2n_reindex(old_data,g::AdaptivityGlue) -> new_data`
+
+  Reindexes a cell-wise array from the old mesh to the new mesh.
+"""
+function o2n_reindex(old_data,g::AdaptivityGlue{GT,Dc}) where {GT,Dc}
+  ncell_to_ocell = g.n2o_faces_map[Dc+1]
+  return _reindex(old_data,ncell_to_ocell)
 end
 
-f2c_reindex(a::CellDatum,g::AdaptivityGlue) = f2c_reindex(CellData.get_data(a),g::AdaptivityGlue)
-c2f_reindex(a::CellDatum,g::AdaptivityGlue) = c2f_reindex(CellData.get_data(a),g::AdaptivityGlue)
+n2o_reindex(a::CellDatum,g::AdaptivityGlue) = n2o_reindex(CellData.get_data(a),g::AdaptivityGlue)
+o2n_reindex(a::CellDatum,g::AdaptivityGlue) = o2n_reindex(CellData.get_data(a),g::AdaptivityGlue)
 
 function _reindex(data,idx::Table)
   m = Reindex(data)
@@ -189,7 +216,7 @@ function get_d_to_fface_to_cface(glue::AdaptivityGlue{<:RefinementGlue},
                                  ftopo::GridTopology{Dc}) where Dc
 
   # Local data for each coarse cell, at the RefinementRule level
-  rrules = Adaptivity.get_old_cell_refinement_rules(glue)
+  rrules = 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)
 

src/Adaptivity/RefinementRules.jl

@@ -3,6 +3,15 @@ abstract type RefinementRuleType end
 struct GenericRefinement <: RefinementRuleType end
 struct WithoutRefinement <: RefinementRuleType end
 
+"""
+  Structure representing the map between a single parent cell and its children.
+
+  Contains: 
+
+  - T :: `RefinementRuleType`, indicating the refinement method.
+  - poly :: `Polytope`, representing the geometry of the parent cell.
+  - ref_grid :: `DiscreteModel` defined on `poly`, giving the parent-to-children cell map. 
+"""
 struct RefinementRule{P,A<:DiscreteModel}
   T         :: RefinementRuleType
   poly      :: P

test/AdaptivityTests/runtests.jl

@@ -2,24 +2,20 @@ module AdaptivityTests
 
 using Test
 
-@testset "AdaptedGeometry" begin
+@testset "Refinement" begin
   include("RefinementRulesTests.jl")
   include("AdaptedGeometryTests.jl")
-end
-
-@testset "Refinement" begin
+  include("FaceLabelingTests.jl")
   include("CartesianRefinementTests.jl")
   include("ComplexChangeDomainTests.jl")
   include("EdgeBasedRefinementTests.jl")
   include("FineToCoarseFieldsTests.jl")
+  include("RefinementRuleBoundaryTests.jl")
+  include("MultifieldRefinementTests.jl")
 end
 
 @testset "CompositeQuadratures" begin
   include("CompositeQuadratureTests.jl")
 end
 
-@testset "MultiFields" begin
-  include("MultifieldRefinementTests.jl")
-end
-
 end # module