Add newest vertex bisection into new adaptivity framework

[aerappa]

Add the capability (limited to TRI in 2D) to refine based on newest vertex (longest edge) bisection. This gives rise to a small number of self-similar triangles so that, if the initial mesh is well conditioned, all resulting meshes will be also. To achieve this I added the standard blue refinement in red/green/blue refinement as well as a "double blue" refinement which is essentially applying blue refinement twice in one iteration. This is necessary when three edges are marked.

The changes are all contained in src/Adaptivity/EdgeBasedRefinement.jl and test/AdaptivityTests/EdgeBasedRefinementTests.jl. The other changes are just removing the old version from src/Geometry.

@JordiManyer @amartinhuertas

[codecov-commenter]

Codecov Report

Merging #901 (78613f6) into master (af46eea) will decrease coverage by 0.02%.
The diff coverage is 88.37%.

❗ Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

@@            Coverage Diff             @@
##           master     #901      +/-   ##
==========================================
- Coverage   88.51%   88.50%   -0.02%     
==========================================
  Files         173      172       -1     
  Lines       20269    20121     -148     
==========================================
- Hits        17942    17808     -134     
+ Misses       2327     2313      -14     

Impacted Files	Coverage Δ
src/Adaptivity/AdaptedDiscreteModels.jl	83.72% <83.33%> (-0.43%)	⬇️
src/Adaptivity/EdgeBasedRefinement.jl	93.84% <88.55%> (-3.55%)	⬇️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[JordiManyer]

Hey @aerappa , thank you for the all the good work!

I've had a look at the PR, and have one main doubt/suggestion:

I don't see why

abstract type RefinementStrategy end

is needed at all. Isn't it redundant with the already existent types? I.e

abstract type AdaptivityMethod end
struct EdgeBasedRefinement <: AdaptivityMethod end

If you want to add some more complexity, I would do it like so:

abstract type AdaptivityMethod end
abstract type EdgeBasedRefinement <: AdaptivityMethod end

struct NVBRefinement{T} <: EdgeBasedRefinement end
struct RedGreenRefinement <: EdgeBasedRefinement end

This gets rid of the new type, while still allowing you to do the same dispatching.

[JordiManyer]

With the above change, some other things also get cleaner, notably

function refine(::EdgeBasedRefinement,model::UnstructuredDiscreteModel{Dc,Dp};cells_to_refine=nothing, should_use_nvb=false) where {Dc,Dp}
  # cells_to_refine can be
  #    a) nothing -> All cells get refined
  #    b) AbstractArray{<:Bool} of size num_cells(model)
  #            -> Only cells such that cells_to_refine[iC] == true get refined
  #    c) AbstractArray{<:Integer}
  #            -> Cells for which gid ∈ cells_to_refine get refined

  # Create new model
  if should_use_nvb
    strategy = NVBStrategy(model)
  else
    strategy = RedGreenStrategy()
  end
  rrules, faces_list = setup_edge_based_rrules(strategy, model.grid_topology,cells_to_refine)
  #rrules, faces_list = setup_edge_based_rrules(strategy, model.grid_topology,cells_to_refine)
  topo   = _refine_unstructured_topology(model.grid_topology,rrules,faces_list)
  reffes = map(p->LagrangianRefFE(Float64,p,1),get_polytopes(topo))
  grid   = UnstructuredGrid(get_vertex_coordinates(topo),get_faces(topo,Dc,0),reffes,get_cell_type(topo),OrientationStyle(topo))
  labels = FaceLabeling(topo)
  ref_model = UnstructuredDiscreteModel(grid,topo,labels)
  ## Create ref glue
  glue = _get_refinement_glue(topo,model.grid_topology,rrules)

  return AdaptedDiscreteModel(ref_model,model,glue)
end

becomes

function refine(method::EdgeBasedRefinement,model::UnstructuredDiscreteModel{Dc,Dp};cells_to_refine=nothing) where {Dc,Dp}
  # cells_to_refine can be
  #    a) nothing -> All cells get refined
  #    b) AbstractArray{<:Bool} of size num_cells(model)
  #            -> Only cells such that cells_to_refine[iC] == true get refined
  #    c) AbstractArray{<:Integer}
  #            -> Cells for which gid ∈ cells_to_refine get refined

  # Create new model
  rrules, faces_list = setup_edge_based_rrules(method, model.grid_topology,cells_to_refine)
  topo   = _refine_unstructured_topology(model.grid_topology,rrules,faces_list)
  reffes = map(p->LagrangianRefFE(Float64,p,1),get_polytopes(topo))
  grid   = UnstructuredGrid(get_vertex_coordinates(topo),get_faces(topo,Dc,0),reffes,get_cell_type(topo),OrientationStyle(topo))
  labels = FaceLabeling(topo)
  ref_model = UnstructuredDiscreteModel(grid,topo,labels)
  ## Create ref glue
  glue = _get_refinement_glue(topo,model.grid_topology,rrules)

  return AdaptedDiscreteModel(ref_model,model,glue)
end

which uses dispatching instead of flags to choose between methods. I find this goes better with Julia's design.

[JordiManyer]

Other than this, there are a couple typos in your docstrings, as well as some leftover debgging commented code we should probably remove.

Could you apply the changes and tag me again so that I can have a final check?
Thanks!

[amartinhuertas]

@aerappa Thanks a lot for your work! Can you please add an entry in the NEWS.md file with a succint description of the PR?

[JordiManyer]

Thanks @aerappa for the changes!

I've been considering how future-proof the current interface is as we progressively add more AdaptivityMethods, and I think the current implementation is not general enough (which is definitely my fault from when I implemented this...).

We currently have

function refine(model::UnstructuredDiscreteModel,::AdaptivityMethod,args...;kwargs...)
  @abstractmethod
end

function refine(model::UnstructuredDiscreteModel,args...;should_use_nvb=false,kwargs...)
  # These types are defined in EdgeBasedRefinement.jl
  if should_use_nvb
    refinement_method = NVBRefinement(model)
  else
    refinement_method = RedGreenRefinement()
  end
  return refine(refinement_method,model,args...;kwargs...)
end

which does not really allow for more methods to be built into it in an organic way.

I think a good solution would be the following

function refine(::AdaptivityMethod,model::UnstructuredDiscreteModel,args...;kwargs...)
  @abstractmethod
end

function refine(model::UnstructuredDiscreteModel,args...;refinement_method=RedGreenRefinement(),kwargs...)
  return refine(refinement_method,model,args...;kwargs...)
end

This has the slight downside of the user having to call NVBRefinement(model) outside the function, but I think we can live with that. How do you see it? I am open to suggestions on how to solve this in an elegant way.

[JordiManyer]

A minor comment: If I understand correctly, NVB only works for triangles, right? We could add a @check that ensure that.
i think the NVBRefinement constructor would be a good place to put it, i.e

function NVBRefinement(model::UnstructuredDiscreteModel{Dc,Dp}) where {Dc, Dp}
  topo = model.grid_topology
  @check all(p->p==TRI,get_cell_polytopes(topo))
  c2e_map     = get_faces(topo,Dc,1)
  e2n_map     = get_faces(topo,1 ,0)
  node_coords = get_node_coordinates(model)
  longest_edge_lids, longest_edge_gids = _get_longest_edge_ids(c2e_map, e2n_map, node_coords)
  NVBRefinement(longest_edge_lids, longest_edge_gids)
end

what do you think?

[aerappa]

Hey @JordiManyer. Good point, having the flag as a keyword argument is a little clunky, but I was doing this exactly because I didn't want to create a new NVBRefinement as the user 😅. It's really not a big deal though since it can even be created directly in the call à la refine(model, refinement_method=NVBRefinement(model)). I'll update this in a new commit.

Good point about the check. That is correct it's limited to TRI.

By the way, what do you think about the terminology newest vertex bisection? Historically people called it that, but it's actually more accurate to say longest edge bisection since the algorithm always bisects the longest edge. Anyway, it's just a detail.

[JordiManyer]

Hey @aerappa .

Ok, then let's leave it as so. We'll call it from the outside.
Concerning the terminology, it's your call. I think we should probably follow the literature, even if it is a little contradictory. But LVB looks ok to me as well.

[aerappa]

Hey @JordiManyer, I came up with (what I hope is) a compromise. The user inputs a string from a set list and there is a lookup internally for which type to use.

# Handle the user's requested choice for refinement
function string_to_refinement(refinement_method::String, model)
  refinement_method == "red_green" && return RedGreenRefinement()
  refinement_method == "nvb" && return NVBRefinement(model)
  error("refinement_method $refinement_method not recognized")
end

function refine(model::UnstructuredDiscreteModel,args...;refinement_method="red_green",kwargs...)
  return refine(string_to_refinement(refinement_method, model),model,args...;kwargs...)
end

I personally find this cleaner because NVBRefinement doesn't have to be exported and the fact that NVBRefinement constructor takes the model is hidden. Also it "should" be easily extensible by adding a new string/type pair. But this should probably be documented somewhere.. what do you think? In the doc for the refine function?

[JordiManyer]

Hey @aerappa , I think that's a good compromise. If you are happy with the PR, I am as well. Once you finalize things, tell Alberto and he'll close the PR.

[aerappa]

Hey @amartinhuertas, from my end everything is good here.

[amartinhuertas]

ok, thanks for your work @aerappa ! Accepting PR.