Newest Vertex Bisection

[aerappa]

This is a feature request for the newest vertex bisection algorithm. For a very nice overview of the algorithm including MATLAB code see this paper by Long Chen.

Based on preliminary discussion, the interface for generating one level of refinement might be something like

model_ref = newest_vertex_bisection(model, marked; hanging_allowed=false)

where model is the current discrete model, marked is a list of cells to refine and hanging_allowed is a flag to specify if conformity must be preserved.

@fverdugo has kindly provided a skeleton version of how this might be achieved:

model = # whatherver discrete model you want
grid = get_grid(model)
node_coords = get_node_coordinates(grid)
cell_node_ids  = get_cell_dof_ids(grid)
reffes = get_reffes(grid)

# Now compute your "refined" quantities
nodal_coords_ref = 
cell_node_ids_ref =

# Build the refined grid
ncells_ref = length(cell_node_ids_ref)
cell_type_ref = fill(1,ncells_ref )
 grid_ref = UnstructuredGrid(nodal_coords_ref,cell_node_ids_ref,reffes,cell_type_ref)

# Build the model
topo_ref = GridTopology(grid_ref)
labels_ref = FaceLabeling(topo_ref)
model_ref = DiscreteModel(grid_ref,topo_ref,labels_ref)

This would ideally fit in as the last step in the standard

SOLVE -> ESTIMATE -> MARK -> REFINE

loop. Thus, FEFunctions would need to be interpolated (hopefully in a rather lightweight way) between the various mesh levels of the refinement process. My initial thought was that this may be achievable using BackgroundTriangulations since the algorithm produces a hierarchical set of Triangulations i.e. Tk​⊂T{k+1} for all k.​

[fverdugo]

Hi @aerappa!

The low-level functionality is clear. Now, we need to think a nice high-level API to express the solve/estimate/mark/refine loop in such a way that we can consider several different refinement strategies. Let me think a bit a bout this.

[fverdugo]

Hi @aerappa

as a first step, we can implement a simpler version of the problem and add more complexity later. I think that it can also be a good exercise to learn more Gridap details.

So for the moment I would do:

1. Just refinement of a grid
2. Just refinement of a model (preserving face tags). Except the part related with the face tags, the rest should be easy once the method in step 1) is available.

# setp 1
function newest_vertex_bisection(grid::Grid,cell_mask::AbstractVector{<:Bool})
  # todo
  ref_grid
end

# step 2
function newest_vertex_bisection(model::DiscreteModel,cell_mask::AbstractVector{<:Bool})
  grid  = get_grid(model)
  ref_grid = newest_vertex_bisection(grid,cell_mask)
  ref_topo = GridTopology(grid)
  labels = get_face_labelling(model)
  ref_labels = # Compute them from the original labels (This is perhaps the most tedious part)
  ref_model = DiscreteModel(ref_grid,ref_topo,ref_labels)
  ref_model
end

At a later stage, we can do the projection onto the new grid etc.

[fverdugo]

And for the moment, I would refine ALL cells so that you don't have hanging nodes.

FESpaces on models with hanging nodes is still not implemented in Gridap, but we are close. We already have the FESpaceWIthLinearConstraints that will to the dirty work.

[aerappa]

Hi @fverdugo I finally got some time to start working on this. For now I am just looking at 2d. I noticed that the way cells are stored according to either

model = some model..
 grid  = get_grid(model)
cell_node_ids = get_cell_node_ids(grid)

or (seemingly) equivalently

top = get_grid_topology(model)
cell_node_ids = get_faces(top, 2, 0)

appear to be just sorted according to the global id of the nodes that define them. It would be more convenient for bisection algorithms if I had a data structure where the cells were sorted in a locally consistent way i.e. clockwise or anti-clockwise. This is useful when constructing the dual graph. I am wondering if such a sorted version of cell_node_ids is possible to extract or if I need to write a conversion layer that uses the underlying geometry to sort?

Edit: I have created an interface to sort each entry of cell_node_ids anti-clockwise. However, do I need to resort this array based on GIDs before creating a new UnstructuredGrid?

[fverdugo]

Edit: I have created an interface to sort each entry of cell_node_ids anti-clockwise. However, do I need to resort this array based on GIDs before creating a new UnstructuredGrid?

Which reference finite elements are using for your FESpace? Lagrangian elements on simpices should work fine with any ordering. So, I think that you don't need to reorder.

BTW,

cell_node_ids = get_cell_node_ids(grid)

and

cell_vertex_ids = get_faces(top, 2, 0)

are not equivalent in general. One gives node ids and the other vertex ids. For linear grids nodes and vertices are equivalent, but for high-order grids they are different.

[fverdugo]

@amartinhuertas, is the orientation based on global nodes ids needed any more after your refactoring of the Raviart-Thomas reference element?

If not, we can use a more standard orientation for simplex meshes, ie same sign of det(Jt) in all cells. I think it would simplify things.

cc @santiagobadia

[aerappa]

@fverdugo It seems the ordering based on GID doesn't matter (or maybe it's getting sorted internally somewhere upon calling UnstructuredGrid). I have a basic version of the algorithm working. See attached image for refinement of all the cells on a 3x3 simplicial mesh. Should I make a PR so that we can discuss the code? For instance, I was actually able to avoid manually re-creating the ref_labels by calling FaceLabeling(ref_topo) but I'm sure this wouldn't work in general. Update: Nonuniform refinement seems ok too with no hanging nodes