From e3da0a4fb4eb5809402f57d52fff80344e035478 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 27 Jun 2024 16:49:14 +1000
Subject: [PATCH 1/8] Minor change in get_o2n_faces_map

---
 src/Adaptivity/AdaptivityGlues.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/Adaptivity/AdaptivityGlues.jl b/src/Adaptivity/AdaptivityGlues.jl
index b94dde743..76015a18b 100644
--- a/src/Adaptivity/AdaptivityGlues.jl
+++ b/src/Adaptivity/AdaptivityGlues.jl
@@ -122,15 +122,15 @@ function get_o2n_faces_map(ncell_to_ocell::Vector{T}) where {T<:Integer}
     iC = ncell_to_ocell[iF]
     ptrs[iC+1] += 1
   end
-  length_to_ptrs!(ptrs)
+  Arrays.length_to_ptrs!(ptrs)
 
-  cnts = fill(0,nC)
   data = fill(zero(T),ptrs[end])
   for iF in 1:nF
     iC = ncell_to_ocell[iF]
-    data[ptrs[iC]+cnts[iC]] = iF
-    cnts[iC] += 1
+    data[ptrs[iC]] = iF
+    ptrs[iC] += 1
   end
+  Arrays.rewind_ptrs!(ptrs)
 
   ocell_to_ncell = Table(data,ptrs)
   return ocell_to_ncell

From c4067b799d61a3b558865b265428b0b928a87313 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Mon, 8 Jul 2024 17:22:54 +1000
Subject: [PATCH 2/8] 3D, Barycentric and Powell Sabin Refinement

---
 src/Adaptivity/AdaptedDiscreteModels.jl |   2 +-
 src/Adaptivity/AdaptivityGlues.jl       |   6 +-
 src/Adaptivity/EdgeBasedRefinement.jl   | 552 +++++++++++++++++-------
 src/Arrays/Tables.jl                    |   6 +
 4 files changed, 416 insertions(+), 150 deletions(-)

diff --git a/src/Adaptivity/AdaptedDiscreteModels.jl b/src/Adaptivity/AdaptedDiscreteModels.jl
index 0b91d51c0..3a640d835 100644
--- a/src/Adaptivity/AdaptedDiscreteModels.jl
+++ b/src/Adaptivity/AdaptedDiscreteModels.jl
@@ -134,7 +134,7 @@ function refine(model::CartesianDiscreteModel{Dc}, cell_partition::Tuple) where
   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)
+  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)
diff --git a/src/Adaptivity/AdaptivityGlues.jl b/src/Adaptivity/AdaptivityGlues.jl
index 76015a18b..bfb434169 100644
--- a/src/Adaptivity/AdaptivityGlues.jl
+++ b/src/Adaptivity/AdaptivityGlues.jl
@@ -261,17 +261,17 @@ end
 
 # FaceLabeling refinement
 
-function _refine_face_labeling(coarse_labeling::FaceLabeling,
+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)
+  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,
+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)
diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl
index 53452b995..83b9f2c7d 100644
--- a/src/Adaptivity/EdgeBasedRefinement.jl
+++ b/src/Adaptivity/EdgeBasedRefinement.jl
@@ -29,57 +29,20 @@ Note on RefinementRules and Orientation of the refined grids:
       X_coarse = Φ∘β(X_fine)
 """
 
-abstract type EdgeBasedRefinement <: AdaptivityMethod end
-
 
-struct NVBRefinement{T} <: EdgeBasedRefinement
-  cell_to_longest_edge_lid::Vector{T}
-  cell_to_longest_edge_gid::Vector{T}
-end
-
-_calculate_edge_length(node_coords, nodes) = norm(node_coords[nodes[1]] - node_coords[nodes[2]])
+"""
+    abstract type EdgeBasedRefinement <: AdaptivityMethod end
 
-# Create Vectors containing the longest edge gids and lids
-function _get_longest_edge_ids(c2e_map, e2n_map, node_coords)
-  nC = length(c2e_map)
-  e2n_map_cache = array_cache(e2n_map)
-  c2e_map_cache = array_cache(c2e_map)
-  longest_edge_gids = Vector{eltype(eltype(c2e_map))}(undef, nC)
-  longest_edge_lids = Vector{eltype(eltype(c2e_map))}(undef, nC)
-  # Loop over cells
-  for c in 1:length(c2e_map)
-    es = getindex!(c2e_map_cache, c2e_map, c)
-    e_length_max = 0.0
-    e_gid_max = -1 # Longest edge index must be found
-    e_lid_max = -1 # Longest edge index must be found
-    for (e_lid, e_gid) in enumerate(es)
-      ns = getindex!(e2n_map_cache, e2n_map, e_gid)
-      e_length = _calculate_edge_length(node_coords, ns)
-      if e_length > e_length_max
-        e_length_max = e_length
-        e_gid_max = e_gid
-        e_lid_max = e_lid
-      end
-    end
-    longest_edge_gids[c] = e_gid_max
-    longest_edge_lids[c] = e_lid_max
-  end
-  return longest_edge_lids, longest_edge_gids
-end
+Abstract type representing edge-based refinement methods. Edge-based refinement methods
+are the ones where all new nodes are created from the faces of the parent mesh.
+This allows the new connectivity to easily be computed from the parent connectivity.
 
-function NVBRefinement(model::DiscreteModel{Dc,Dp}) where {Dc, Dp}
-  topo = model.grid_topology
-  @check all(p->p==TRI,get_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
+Any edge-based refinement method has to implement the following API: 
 
-NVBRefinement(model::AdaptedDiscreteModel) = NVBRefinement(model.model)
+- [`setup_edge_based_rrules`](@ref)
 
-struct RedGreenRefinement <: EdgeBasedRefinement end
+"""
+abstract type EdgeBasedRefinement <: AdaptivityMethod end
 
 function refine(method::EdgeBasedRefinement,model::UnstructuredDiscreteModel{Dc,Dp};cells_to_refine=nothing) where {Dc,Dp}
   # cells_to_refine can be
@@ -92,37 +55,38 @@ function refine(method::EdgeBasedRefinement,model::UnstructuredDiscreteModel{Dc,
   coarse_labels = get_face_labeling(model)
   # Create new model
   rrules, faces_list = setup_edge_based_rrules(method, model.grid_topology,cells_to_refine)
-  topo   = _refine_unstructured_topology(ctopo,rrules,faces_list)
+  topo   = refine_edge_based_topology(ctopo,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))
-  glue = _get_refinement_glue(topo,model.grid_topology,rrules)
-  labels = _refine_face_labeling(coarse_labels,glue,model.grid_topology,topo)
+  grid   = UnstructuredGrid(
+    get_vertex_coordinates(topo),
+    get_faces(topo,Dc,0),reffes,
+    get_cell_type(topo),
+    OrientationStyle(topo)
+  )
+  glue = get_edge_based_refinement_glue(topo,model.grid_topology,rrules)
+  labels = refine_face_labeling(coarse_labels,glue,model.grid_topology,topo)
   ref_model = UnstructuredDiscreteModel(grid,topo,labels)
   return AdaptedDiscreteModel(ref_model,model,glue)
 end
 
-function _refine_unstructured_topology(topo::UnstructuredGridTopology{Dc},
-                                      rrules::AbstractVector{<:RefinementRule},
-                                      faces_list::Tuple) where {Dc}
+function refine_edge_based_topology(
+  topo::UnstructuredGridTopology{Dc},
+  rrules::AbstractVector{<:RefinementRule},
+  faces_list::Tuple
+) where {Dc}
   coords_new  = get_new_coordinates_from_faces(topo,faces_list)
   c2n_map_new = get_refined_cell_to_vertex_map(topo,rrules,faces_list)
   polys_new, cell_type_new = _get_cell_polytopes(rrules)
-
-  # We can guarantee the new topology is oriented if
-  #   1 - the old topology was oriented
-  #   2 - we have a single type of polytope (i.e new topo is not mixed)
-  #orientation = NonOriented()
   orientation = NonOriented()
 
   return UnstructuredGridTopology(coords_new,c2n_map_new,cell_type_new,polys_new,orientation)
 end
 
-function _get_refinement_glue(ftopo ::UnstructuredGridTopology{Dc},
-                              ctopo ::UnstructuredGridTopology{Dc},
-                              rrules::AbstractVector{<:RefinementRule}) where {Dc}
+function get_edge_based_refinement_glue(
+  ftopo ::UnstructuredGridTopology{Dc},
+  ctopo ::UnstructuredGridTopology{Dc},
+  rrules::AbstractVector{<:RefinementRule}
+) where {Dc}
   nC_old = num_faces(ctopo,Dc)
   nC_new = num_faces(ftopo,Dc)
 
@@ -143,27 +107,173 @@ function _get_refinement_glue(ftopo ::UnstructuredGridTopology{Dc},
 end
 
 """
-Given an UnstructuredTopology and a list of cells_to_refine, provides an edge-based coloring
-for the topology:
+    setup_edge_based_rrules(method::EdgeBasedRefinement,topo::UnstructuredGridTopology,cells_to_refine)
 
-  - If a cell is completely refined, it is colored RED
-  - If a cell touches more than one refined (RED) cell, it becomes RED.
-  - If a cell touches a single refined (RED) cell, it is colored GREEN.
+Given an UnstructuredTopology and a list of cells to refine, returns a vector
+of refinement rules and a list of the faces (called refined faces) where new 
+vertices will be created.
+"""
+function setup_edge_based_rrules(::EdgeBasedRefinement,topo::UnstructuredGridTopology{Dc},cells_to_refine) where Dc
+  @abstractmethod
+end
+
+function setup_edge_based_rrules(method::EdgeBasedRefinement,topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Bool}) where Dc
+  return setup_edge_based_rrules(method,topo,findall(cells_to_refine))
+end
 
-The method returns a vector of Red/Green refinement rules, as well the list of
-vertices, edges and cells which correspond to new vertices in the refined topology.
 """
-function setup_edge_based_rrules(::EdgeBasedRefinement, topo::UnstructuredGridTopology{Dc},::Nothing) where Dc
-  rrules     = lazy_map(RedRefinementRule,CompressedArray(topo.polytopes,topo.cell_type))
-  faces_list = _redgreen_refined_faces_list(topo,rrules,[true])
-  return rrules, faces_list
+    get_new_coordinates_from_faces(p::Union{Polytope{D},GridTopology{D}},faces_list::Tuple)
+
+Given a Polytope or GridTopology and a list of refined d-faces, returns the new coordinates
+of the vertices created on those faces: 
+  - If the face is a 0-face (node), the new vertex is the node itself.
+  - For d > 0, the new vertex is the barycenter of the face.
+The new vertices are ordered by parent dimension and face id (in that order).
+"""
+function get_new_coordinates_from_faces(p::Union{Polytope{D},GridTopology{D}},faces_list::Tuple) where {D}
+  @check length(faces_list) == D+1
+
+  nN_new     = sum(x->length(x),faces_list)
+  coords_old = get_vertex_coordinates(p)
+  coords_new = Vector{eltype(coords_old)}(undef,nN_new)
+
+  n = 1
+  for (d,dfaces) in enumerate(faces_list)
+    if length(dfaces) > 0
+      nf = length(dfaces)
+      d2n_map = get_faces(p,d-1,0)
+      coords_new[n:n+nf-1] .= map(f -> sum(coords_old[d2n_map[f]])/length(d2n_map[f]), dfaces)
+      n += nf
+    end
+  end
+
+  return coords_new
 end
 
-function setup_edge_based_rrules(::EdgeBasedRefinement,topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Bool}) where Dc
-  return setup_edge_based_rrules(topo,findall(cells_to_refine))
+"""
+    get_refined_cell_to_vertex_map(
+      topo::UnstructuredGridTopology{Dc},
+      rrules::AbstractVector{<:RefinementRule},
+      faces_list::Tuple
+    )
+
+Given an UnstructuredGridTopology, a vector of RefinementRules and a list of refined faces,
+returns the cell-to-vertex connectivity of the refined topology.
+"""
+function get_refined_cell_to_vertex_map(
+  topo::UnstructuredGridTopology{Dc},
+  rrules::AbstractVector{<:RefinementRule},
+  faces_list::Tuple
+) where Dc
+  @check length(faces_list) == Dc+1
+  d_to_nfaces = map(Df -> num_faces(topo,Df), 0:Dc)
+  d_to_cell_to_dface = map(Df -> get_faces(topo,Dc,Df), 0:Dc-1)
+
+  # Offsets for new vertices
+  d_to_offset = [1,map(length,faces_list)...]
+  length_to_ptrs!(d_to_offset)
+  d_to_offset = d_to_offset[1:end-1]
+
+  # Maps from refined faces (i.e new vertices) to old faces
+  d_to_faces_reindexing = map(d_to_nfaces,d_to_offset,faces_list) do nfaces,offset,ref_faces
+    faces_reindexing = find_inverse_index_map(ref_faces,nfaces)
+    return map(f -> iszero(f) ? 0 : f + offset, faces_reindexing)
+  end
+
+  # For a given parent cell, return the new node gids (corresponding to the refined faces)
+  function cell_to_reindexed_faces(::Val{2},cell)
+    N = view(d_to_cell_to_dface[1],cell)
+    E = view(d_to_faces_reindexing[2],view(d_to_cell_to_dface[2],cell))
+    C = [d_to_faces_reindexing[3][cell]]
+    return (N,E,C)
+  end
+  function cell_to_reindexed_faces(::Val{3},cell)
+    N = view(d_to_cell_to_dface[1],cell)
+    E = view(d_to_faces_reindexing[2],view(d_to_cell_to_dface[2],cell))
+    F = view(d_to_faces_reindexing[3],view(d_to_cell_to_dface[3],cell))
+    C = [d_to_faces_reindexing[4][cell]]
+    return (N,E,F,C)
+  end
+
+  # Allocate ptr and data arrays for new connectivity
+  nC_new    = sum(rr -> num_subcells(rr), rrules)
+  nData_new = sum(rr -> sum(length,rr.ref_grid.grid.cell_node_ids), rrules)
+  ptrs_new  = Vector{Int}(undef,nC_new+1)
+  data_new  = Vector{Int}(undef,nData_new)
+
+  k = 1
+  ptrs_new[1] = 1
+  for iC = 1:nC_old
+    rr = rrules[iC]
+
+    # New vertex ids from old face ids
+    faces = cell_to_reindexed_faces(Val(Dc),iC)
+    sub_conn = get_relabeled_connectivity(rr,faces)
+
+    nChild = length(sub_conn)
+    ptrs_new[k:k+nChild] .= sub_conn.ptrs .+ (ptrs_new[k] - 1)
+    data_new[ptrs_new[k]:ptrs_new[k+nChild]-1] .= sub_conn.data
+    k = k+nChild
+  end
+
+  return Table(data_new,ptrs_new)
 end
 
-_shift_to_first(v::AbstractVector{T}, i::T) where {T<:Integer} = circshift(v, -(i - 1))
+############################################################################################
+
+# NVB Refinement (Nearest Vertex Bissection)
+
+"""
+    struct NVBRefinement{T} <: EdgeBasedRefinement
+
+A refinement method for TRIs, where a selected cell is refined by bisecting the longest edge.
+"""
+struct NVBRefinement{T} <: EdgeBasedRefinement
+  cell_to_longest_edge_lid::Vector{T}
+  cell_to_longest_edge_gid::Vector{T}
+end
+
+_calculate_edge_length(node_coords, nodes) = norm(node_coords[nodes[1]] - node_coords[nodes[2]])
+
+# Create Vectors containing the longest edge gids and lids
+function _get_longest_edge_ids(c2e_map, e2n_map, node_coords)
+  nC = length(c2e_map)
+  e2n_map_cache = array_cache(e2n_map)
+  c2e_map_cache = array_cache(c2e_map)
+  longest_edge_gids = Vector{eltype(eltype(c2e_map))}(undef, nC)
+  longest_edge_lids = Vector{eltype(eltype(c2e_map))}(undef, nC)
+  # Loop over cells
+  for c in 1:length(c2e_map)
+    es = getindex!(c2e_map_cache, c2e_map, c)
+    e_length_max = 0.0
+    e_gid_max = -1 # Longest edge index must be found
+    e_lid_max = -1 # Longest edge index must be found
+    for (e_lid, e_gid) in enumerate(es)
+      ns = getindex!(e2n_map_cache, e2n_map, e_gid)
+      e_length = _calculate_edge_length(node_coords, ns)
+      if e_length > e_length_max
+        e_length_max = e_length
+        e_gid_max = e_gid
+        e_lid_max = e_lid
+      end
+    end
+    longest_edge_gids[c] = e_gid_max
+    longest_edge_lids[c] = e_lid_max
+  end
+  return longest_edge_lids, longest_edge_gids
+end
+
+function NVBRefinement(model::DiscreteModel{Dc,Dp}) where {Dc, Dp}
+  topo = model.grid_topology
+  @check all(p->p==TRI,get_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
+
+NVBRefinement(model::AdaptedDiscreteModel) = NVBRefinement(model.model)
 
 function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}) where Dc
   nE = num_faces(topo,1)
@@ -249,10 +359,47 @@ function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTo
   end
   rrules = CompressedArray(color_rrules,cell_color)
   face_list = _redgreen_refined_faces_list(topo, rrules, is_refined)
+  
   return rrules, face_list
 end
 
-function setup_edge_based_rrules(::RedGreenRefinement, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}) where Dc
+############################################################################################
+
+# RedGreen Refinement
+
+"""
+    struct RedGreenRefinement <: EdgeBasedRefinement
+
+Red-Green conformal refinement method. The refinement follows the following rules:
+
+  - If a cell is completely refined, it is colored RED
+  - If a cell touches more than one refined (RED) cell, it becomes RED.
+  - If a cell touches a single refined (RED) cell, it is colored GREEN.
+
+The resulting mesh is always conformal.
+"""
+struct RedGreenRefinement <: EdgeBasedRefinement end
+
+function setup_edge_based_rrules(::RedGreenRefinement, topo::UnstructuredGridTopology{Dc},::Nothing) where Dc
+  rrules     = lazy_map(RedRefinementRule,CompressedArray(topo.polytopes,topo.cell_type))
+
+  ref_nodes = 1:num_faces(topo,0)
+  ref_edges = 1:num_faces(topo,1)
+  ref_cells = findall(_has_interior_point,rrules)
+
+  if Dc == 2
+    faces_list = (ref_nodes,ref_edges,ref_cells)
+  else
+    faces_list = (ref_nodes,ref_edges,Int32[],ref_cells)
+  end
+
+  return rrules, faces_list
+end
+
+function setup_edge_based_rrules(
+  ::RedGreenRefinement, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}
+) where Dc
+  @check Dc == 2
   nC = num_cells(topo)
   nE = num_faces(topo,1)
   c2e_map     = get_faces(topo,Dc,1)
@@ -328,79 +475,53 @@ function _redgreen_refined_faces_list(topo::UnstructuredGridTopology{2},rrules,e
   return (ref_nodes,ref_edges,ref_cells)
 end
 
-function get_new_coordinates_from_faces(p::Union{Polytope{D},GridTopology{D}},faces_list::Tuple) where {D}
-  @check length(faces_list) == D+1
+############################################################################################
 
-  nN_new     = sum(x->length(x),faces_list)
-  coords_old = get_vertex_coordinates(p)
-  coords_new = Vector{eltype(coords_old)}(undef,nN_new)
+# Barycentric refinement
 
-  n = 1
-  for (d,dfaces) in enumerate(faces_list)
-    if length(dfaces) > 0
-      nf = length(dfaces)
-      d2n_map = get_faces(p,d-1,0)
-      coords_new[n:n+nf-1] .= map(f -> sum(coords_old[d2n_map[f]])/length(d2n_map[f]), dfaces)
-      n += nf
-    end
-  end
-
-  return coords_new
-end
+"""
+    struct BarycentricRefinementMethod <: EdgeBasedRefinement end
 
-function get_refined_cell_to_vertex_map(topo::UnstructuredGridTopology{2},
-                                        rrules::AbstractVector{<:RefinementRule},
-                                        faces_list::Tuple)
-  @notimplementedif !all(map(p -> p ∈ [QUAD,TRI],topo.polytopes))
-  nN_old,nE_old,nC_old  = num_faces(topo,0),num_faces(topo,1),num_faces(topo,2)
-  c2n_map = get_faces(topo,2,0)
-  c2e_map = get_faces(topo,2,1)
-  ref_edges = faces_list[2]
+Barycentric refinement method. Each selected cell is refined by adding a new vertex at the
+barycenter of the cell. Then the cell is split into subcells by connecting the new vertex
+to the vertices of the parent cell. Always results in a conformal mesh.
+"""
+struct BarycentricRefinementMethod <: EdgeBasedRefinement end
 
-  # Allocate map ptr and data arrays
-  nC_new    = sum(rr -> num_subcells(rr), rrules)
-  nData_new = sum(rr -> sum(ids->length(ids),rr.ref_grid.grid.cell_node_ids), rrules)
+function setup_edge_based_rrules(
+  method::BarycentricRefinementMethod, topo::UnstructuredGridTopology,::Nothing
+)
+  setup_edge_based_rrules(method,topo,collect(1:num_faces(topo,Dc)))
+end
 
-  ptrs_new  = Vector{Int}(undef,nC_new+1)
-  data_new  = Vector{Int}(undef,nData_new)
+function setup_edge_based_rrules(
+  ::BarycentricRefinementMethod, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}
+) where Dc
+  rrules = lazy_map(BarycentricRefinementRule,CompressedArray(topo.polytopes,topo.cell_type))
 
-  cell_is_refined = lazy_map(rr->isa(RefinementRuleType(rr),RedRefinement),rrules)
-  if all(cell_is_refined)
-    edges_reindexing = 1:nE_old
+  ref_nodes = 1:num_faces(topo,0)
+  if Dc == 2
+    faces_list = (ref_nodes,Int32[],cells_to_refine)
   else
-    edges_reindexing = lazy_map(Reindex(find_inverse_index_map(ref_edges,nE_old)),1:nE_old)
-  end
-
-  k = 1
-  ptrs_new[1] = 1
-  C = nN_old + length(ref_edges) + 1
-  for iC = 1:nC_old
-    rr = rrules[iC]
-
-    # New Node gids from old N,E,C lids
-    N = c2n_map[iC]
-    E = edges_reindexing[c2e_map[iC]] .+ nN_old
-    sub_conn = get_relabeled_connectivity(rr,(N,E,[C]))
-
-    nChild = length(sub_conn)
-    ptrs_new[k:k+nChild] .=  sub_conn.ptrs .+ (ptrs_new[k] - 1)
-    data_new[ptrs_new[k]:ptrs_new[k+nChild]-1] .= sub_conn.data
-    k = k+nChild
-
-    _has_interior_point(rr) && (C += 1)
+    faces_list = (ref_nodes,Int32[],Int32[],cells_to_refine)
   end
 
-  return Table(data_new,ptrs_new)
+  return rrules, faces_list
 end
 
-###############################################################################
+############################################################################################
+
 # EdgeBasedRefinementRules
 
 abstract type EdgeBasedRefinementRule <: RefinementRuleType end
+
+# List of available RefinementRules:
 struct RedRefinement <: EdgeBasedRefinementRule end
+struct GreenRefinement{P} <: EdgeBasedRefinementRule end
 struct BlueRefinement{P, Q} <: EdgeBasedRefinementRule end
 struct BlueDoubleRefinement{P} <: EdgeBasedRefinementRule end
-struct GreenRefinement{P} <: EdgeBasedRefinementRule end
+struct BarycentricRefinement <: EdgeBasedRefinementRule end
+struct PowellSabinRefinement <: EdgeBasedRefinementRule end
 
 _has_interior_point(rr::RefinementRule) = _has_interior_point(rr,RefinementRuleType(rr))
 _has_interior_point(rr::RefinementRule,::RefinementRuleType) = false
@@ -418,7 +539,7 @@ Edge-based RefinementRule where a new vertex is added to
 each edge of the original Polytope.
 """
 function RedRefinementRule(p::Polytope)
-  @notimplementedif (p ∉ [TRI,QUAD,HEX])
+  @notimplementedif (p ∉ [TRI,TET,QUAD,HEX])
 
   faces_list = _get_red_refined_faces_list(p)
   coords = get_new_coordinates_from_faces(p,faces_list)
@@ -433,6 +554,8 @@ end
 function _get_red_refined_faces_list(p::Polytope)
   if p == TRI
     return (Int32[1,2,3],Int32[1,2,3],Int32[])
+  elseif p == TET
+    return (Int32[1,2,3,4],Int32[1,2,3,4,5,6],Int32[],Int32[])
   elseif p == QUAD
     return (Int32[1,2,3,4],Int32[1,2,3,4],Int32[1])
   elseif p == HEX
@@ -451,6 +574,19 @@ function _get_red_refined_connectivity(p::Polytope)
                       4,5,6]
     conn_ptrs = Int32[1,4,7,10,13]
     return polys, cell_type, Table(conn_data,conn_ptrs)
+  elseif p == TET
+    polys     = [TET]
+    cell_type = Int32[1,1,1,1,1,1,1,1]
+    conn_data = Int32[1,5,6,8,
+                      5,2,7,9,
+                      6,7,3,10,
+                      8,9,10,4,
+                      5,7,10,9,
+                      8,5,10,9,
+                      5,10,7,6,
+                      8,10,5,6]
+    conn_ptrs = Int32[1,5,9,13,17,21,25,29,33]
+    return polys, cell_type, Table(conn_data,conn_ptrs)
   elseif p == QUAD
     polys     = [QUAD]
     cell_type = Int32[1,1,1,1]
@@ -479,10 +615,7 @@ end
 
 function _has_interior_point(rr::RefinementRule,::RedRefinement)
   p = get_polytope(rr)
-  if p ∈ [QUAD,HEX]
-    return true
-  end
-  return false
+  return !is_simplex(p)
 end
 
 # [Face dimension][Coarse Face id] -> [Fine faces]
@@ -731,6 +864,114 @@ function _get_blue_double_refined_connectivity(p::Polytope{2}, long_ref_edge)
   @notimplemented
 end
 
+"""
+RefinementRule representing barycentric refinement of a Polytope. A single
+vertex is added in the center, then joined to the vertices of the original
+Polytope.
+"""
+function BarycentricRefinementRule(p::Polytope)
+  @notimplementedif (p ∉ [TRI,TET])
+
+  faces_list = _get_barycentric_refined_faces_list(p)
+  coords = get_new_coordinates_from_faces(p,faces_list)
+
+  polys, cell_types, conn = _get_barycentric_refined_connectivity(p)
+  reffes = map(x->LagrangianRefFE(Float64,x,1),polys)
+
+  ref_grid = UnstructuredDiscreteModel(UnstructuredGrid(coords,conn,reffes,cell_types))
+  return RefinementRule(BarycentricRefinement(),p,ref_grid)
+end
+
+function _get_barycentric_refined_faces_list(p::Polytope)
+  if p == TRI
+    return (Int32[1,2,3],Int32[],Int32[1])
+  elseif p == TET
+    return (Int32[1,2,3,4],Int32[],Int32[],Int32[1])
+  end
+  @notimplemented
+end
+
+function _get_barycentric_refined_connectivity(p::Polytope)
+  if p == TRI
+    polys     = [TRI]
+    cell_type = [1, 1, 1]
+    conn_data = [1, 2, 4,
+                 2, 3, 4,
+                 3, 1, 4]
+    conn_ptrs = [1, 4, 7]
+    return polys, cell_type, Table(conn_data,conn_ptrs)
+  elseif p == TET
+    polys     = [TET]
+    cell_type = [1, 1, 1, 1]
+    conn_data = [1, 2, 3, 5,
+                 2, 3, 4, 5,
+                 3, 1, 4, 5,
+                 1, 2, 4, 5]
+    conn_ptrs = [1, 5, 9, 13, 17]
+    return polys, cell_type, Table(conn_data,conn_ptrs)
+  end
+  @notimplemented
+end
+
+"""
+RefinementRule representing the Powell-Sabin split of TRI, and it's extension to
+3-demensional TETs (Worsey-Farin split).
+"""
+function PowellSabinRefinementRule(p::Polytope)
+  @notimplementedif (p ∉ [TRI,TET])
+
+  faces_list = _get_powellsabin_refined_faces_list(p)
+  coords = get_new_coordinates_from_faces(p,faces_list)
+
+  polys, cell_types, conn = _get_powellsabin_refined_connectivity(p)
+  reffes = map(x->LagrangianRefFE(Float64,x,1),polys)
+
+  ref_grid = UnstructuredDiscreteModel(UnstructuredGrid(coords,conn,reffes,cell_types))
+  return RefinementRule(PowellSabinRefinement(),p,ref_grid)
+end
+
+function _get_powellsabin_refined_faces_list(p::Polytope)
+  if p == TRI
+    return (Int32[1,2,3],Int32[1,2,3],Int32[1])
+  elseif p == TET
+    return (Int32[1,2,3,4],Int32[1,2,3,4,5,6],Int32[1,2,3,4],Int32[1])
+  end
+  @notimplemented
+end
+
+function _get_powellsabin_refined_connectivity(p::Polytope)
+  # A) For TRIs, we add a new vertices 
+  #       - in the center of the triangle
+  #       - in the middle of each edge
+  # and then connect the center vertex to all other vertices (nodes + new edge vertices). 
+  # This yields 6 new triangles.
+  # B) For TETs, we do do the same process for each face (which is a TRI), and then connect 
+  #    every new triangle to the center of the TET. This yields 24=n_TRI*n_FACE new TETs.
+
+  n_TRI  = 6 # Number of children in a refined TRI
+  n_FACE = 4 # Number of faces in a TET
+  conn_data_TRI = [1, 4, 7,
+                  4, 2, 7,
+                  2, 6, 7,
+                  6, 3, 7,
+                  3, 5, 7,
+                  5, 1, 7]
+  if p == TRI
+    polys = [TRI]
+    cell_type = fill(1, n_TRI)
+    conn_ptrs = collect(1:3:n_TRI*3+1)
+    return polys, cell_type, Table(conn_data_TRI,conn_ptrs)
+  elseif p == TET
+    polys     = [TET]
+    cell_type = fill(1, n_TRI*n_FACE)
+    face_conn = [[1,2,3,5,6,7,11],[1,2,4,5,8,9,12],[1,3,4,6,8,10,13],[2,3,4,7,9,10,14]]
+    conn_data = vcat([[lazy_map(Reindex(face_conn[f]),conn_data_2d[1+(i-1)*3:i*3])...,15] for i in 1:n_TRI for f in 1:n_FACE]...)
+    conn_ptrs = collect(1:4:n_TRI*n_FACE*4+1)
+    return polys, cell_type, Table(conn_data,conn_ptrs)
+  end
+  @notimplemented
+end
+
 """
 Provided the gids of the coarse cell faces as a Tuple(vertices,edges,cell) and the
 edge-based RefinementRule we want to apply, returns the connectivity of all the refined
@@ -742,35 +983,39 @@ end
 
 get_relabeled_connectivity(::RefinementRuleType,::RefinementRule,faces_gids) = @notimplemented
 
-function get_relabeled_connectivity(::WithoutRefinement,rr::RefinementRule{P},faces_gids) where P<:Polytope{2}
+function get_relabeled_connectivity(::WithoutRefinement,rr::RefinementRule{<:Polytope{2}},faces_gids)
   conn = rr.ref_grid.grid.cell_node_ids
   gids = faces_gids[1]
   new_data = lazy_map(Reindex(gids),conn.data)
   return Table(new_data,conn.ptrs)
 end
 
-function get_relabeled_connectivity(::RedRefinement,rr::RefinementRule{P},faces_gids) where P<:Polytope{2}
+function get_relabeled_connectivity(::RedRefinement,rr::RefinementRule{P},faces_gids) where P
   conn = rr.ref_grid.grid.cell_node_ids
-  gids = [faces_gids[1]...,faces_gids[2]...,faces_gids[3]...]
+  if is_simplex(P) # TRI, TET only adds vertices to edges
+    gids = vcat(faces_gids[1]...,faces_gids[2]...)
+  else # QUAD, HEX adds vertices to edges, faces and interior
+    gids = vcat(faces_gids...)    
+  end
   new_data = lazy_map(Reindex(gids),conn.data)
   return Table(new_data,conn.ptrs)
 end
 
-function get_relabeled_connectivity(::GreenRefinement{N},rr::RefinementRule{P},faces_gids) where {N,P<:Polytope{2}}
+function get_relabeled_connectivity(::GreenRefinement{N},rr::RefinementRule{<:Polytope{2}},faces_gids) where N
   conn = rr.ref_grid.grid.cell_node_ids
   gids = [faces_gids[1]...,faces_gids[2][N]]
   new_data = lazy_map(Reindex(gids),conn.data)
   return Table(new_data,conn.ptrs)
 end
 
-function get_relabeled_connectivity(::BlueRefinement{N, M},rr::RefinementRule{P},faces_gids) where {N, M, P<:Polytope{2}}
+function get_relabeled_connectivity(::BlueRefinement{N, M},rr::RefinementRule{<:Polytope{2}},faces_gids) where {N, M}
   conn = rr.ref_grid.grid.cell_node_ids
   gids = [faces_gids[1]...,faces_gids[2][N], faces_gids[2][M]]
   new_data = lazy_map(Reindex(gids),conn.data)
   return Table(new_data,conn.ptrs)
 end
 
-function get_relabeled_connectivity(::BlueDoubleRefinement{N},rr::RefinementRule{P},faces_gids) where {N, P<:Polytope{2}}
+function get_relabeled_connectivity(::BlueDoubleRefinement{N},rr::RefinementRule{<:Polytope{2}},faces_gids) where {N}
   conn = rr.ref_grid.grid.cell_node_ids
   # Sort needed for non-longest edge gids
   other_ids = setdiff(collect(1:3), N) |> sort
@@ -779,6 +1024,21 @@ function get_relabeled_connectivity(::BlueDoubleRefinement{N},rr::RefinementRule
   return Table(new_data,conn.ptrs)
 end
 
+function get_relabeled_connectivity(::BarycentricRefinement,rr::RefinementRule{<:Polytope{D}},faces_gids) where D
+  conn = rr.ref_grid.grid.cell_node_ids
+  gids = [faces_gids[1]...,faces_gids[D+1]...] # Polytope nodes + Barycenter
+  new_data = lazy_map(Reindex(gids),conn.data)
+  return Table(new_data,conn.ptrs)
+end
+
+function get_relabeled_connectivity(::PowellSabinRefinement,rr::RefinementRule{P},faces_gids) where P
+  @assert is_simplex(P)
+  conn = rr.ref_grid.grid.cell_node_ids
+  gids = vcat(faces_gids...)
+  new_data = lazy_map(Reindex(gids),conn.data)
+  return Table(new_data,conn.ptrs)
+end
+
 function counts_to_ptrs(counts::Vector{<:Integer})
   n = length(counts)
   ptrs = Vector{Int32}(undef,n+1)
diff --git a/src/Arrays/Tables.jl b/src/Arrays/Tables.jl
index 6965635ce..70ef1991e 100644
--- a/src/Arrays/Tables.jl
+++ b/src/Arrays/Tables.jl
@@ -42,6 +42,12 @@ function Base.convert(::Type{Table{T,Vd,Vp}},table::Table{T,Vd,Vp}) where {T,Vd,
   table
 end
 
+function Base.view(a::Table,i::Integer)
+  pini = a.ptrs[i]
+  pend = a.ptrs[i+1]-1
+  return view(a.data,pini:pend)
+end
+
 """
 """
 function identity_table(::Type{T},::Type{P},l::Integer) where {T,P}

From c1d16ac629d55427d6849fdd678ad8194e18e17c Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Mon, 8 Jul 2024 21:55:23 +1000
Subject: [PATCH 3/8] Minor fixes

---
 src/Adaptivity/EdgeBasedRefinement.jl | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl
index 83b9f2c7d..d911663e4 100644
--- a/src/Adaptivity/EdgeBasedRefinement.jl
+++ b/src/Adaptivity/EdgeBasedRefinement.jl
@@ -177,7 +177,7 @@ function get_refined_cell_to_vertex_map(
   # Maps from refined faces (i.e new vertices) to old faces
   d_to_faces_reindexing = map(d_to_nfaces,d_to_offset,faces_list) do nfaces,offset,ref_faces
     faces_reindexing = find_inverse_index_map(ref_faces,nfaces)
-    return map(f -> iszero(f) ? 0 : f + offset, faces_reindexing)
+    return map(f -> iszero(f) ? 0 : f + offset - 1, faces_reindexing)
   end
 
   # For a given parent cell, return the new node gids (corresponding to the refined faces)
@@ -203,7 +203,7 @@ function get_refined_cell_to_vertex_map(
 
   k = 1
   ptrs_new[1] = 1
-  for iC = 1:nC_old
+  for iC = 1:d_to_nfaces[Dc+1]
     rr = rrules[iC]
 
     # New vertex ids from old face ids
@@ -275,6 +275,10 @@ end
 
 NVBRefinement(model::AdaptedDiscreteModel) = NVBRefinement(model.model)
 
+function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTopology{Dc},::Nothing) where Dc
+  setup_edge_based_rrules(method, topo, collect(1:num_faces(topo,Dc)))
+end
+
 function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}) where Dc
   nE = num_faces(topo,1)
   c2e_map       = get_faces(topo,Dc,1)
@@ -381,7 +385,7 @@ The resulting mesh is always conformal.
 struct RedGreenRefinement <: EdgeBasedRefinement end
 
 function setup_edge_based_rrules(::RedGreenRefinement, topo::UnstructuredGridTopology{Dc},::Nothing) where Dc
-  rrules     = lazy_map(RedRefinementRule,CompressedArray(topo.polytopes,topo.cell_type))
+  rrules = lazy_map(RedRefinementRule,CompressedArray(topo.polytopes,topo.cell_type))
 
   ref_nodes = 1:num_faces(topo,0)
   ref_edges = 1:num_faces(topo,1)
@@ -990,9 +994,9 @@ function get_relabeled_connectivity(::WithoutRefinement,rr::RefinementRule{<:Pol
   return Table(new_data,conn.ptrs)
 end
 
-function get_relabeled_connectivity(::RedRefinement,rr::RefinementRule{P},faces_gids) where P
+function get_relabeled_connectivity(::RedRefinement,rr::RefinementRule,faces_gids)
   conn = rr.ref_grid.grid.cell_node_ids
-  if is_simplex(P) # TRI, TET only adds vertices to edges
+  if is_simplex(get_polytope(rr)) # TRI, TET only adds vertices to edges
     gids = vcat(faces_gids[1]...,faces_gids[2]...)
   else # QUAD, HEX adds vertices to edges, faces and interior
     gids = vcat(faces_gids...)    
@@ -1031,8 +1035,8 @@ function get_relabeled_connectivity(::BarycentricRefinement,rr::RefinementRule{<
   return Table(new_data,conn.ptrs)
 end
 
-function get_relabeled_connectivity(::PowellSabinRefinement,rr::RefinementRule{P},faces_gids) where P
-  @assert is_simplex(P)
+function get_relabeled_connectivity(::PowellSabinRefinement,rr::RefinementRule,faces_gids)
+  @assert is_simplex(get_polytope(rr))
   conn = rr.ref_grid.grid.cell_node_ids
   gids = vcat(faces_gids...)
   new_data = lazy_map(Reindex(gids),conn.data)

From 3385ee265f2306a774d3b360256ad687e230d9fd Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Mon, 8 Jul 2024 22:44:35 +1000
Subject: [PATCH 4/8] Added tests

---
 src/Adaptivity/AdaptedDiscreteModels.jl       |  1 +
 src/Adaptivity/EdgeBasedRefinement.jl         | 34 +++++---
 .../EdgeBasedRefinementTests.jl               | 86 +++++++++++++++----
 3 files changed, 93 insertions(+), 28 deletions(-)

diff --git a/src/Adaptivity/AdaptedDiscreteModels.jl b/src/Adaptivity/AdaptedDiscreteModels.jl
index 3a640d835..9c23d81b2 100644
--- a/src/Adaptivity/AdaptedDiscreteModels.jl
+++ b/src/Adaptivity/AdaptedDiscreteModels.jl
@@ -101,6 +101,7 @@ end
 function string_to_refinement(refinement_method::String, model)
   refinement_method == "red_green" && return RedGreenRefinement()
   refinement_method == "nvb" && return NVBRefinement(model)
+  refinement_method == "barycentric" && return BarycentricRefinement()
   error("refinement_method $refinement_method not recognized")
 end
 
diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl
index d911663e4..a60f33b91 100644
--- a/src/Adaptivity/EdgeBasedRefinement.jl
+++ b/src/Adaptivity/EdgeBasedRefinement.jl
@@ -195,6 +195,10 @@ function get_refined_cell_to_vertex_map(
     return (N,E,F,C)
   end
 
+  display(faces_list)
+  display(d_to_faces_reindexing)
+  display(d_to_offset)
+
   # Allocate ptr and data arrays for new connectivity
   nC_new    = sum(rr -> num_subcells(rr), rrules)
   nData_new = sum(rr -> sum(length,rr.ref_grid.grid.cell_node_ids), rrules)
@@ -394,7 +398,10 @@ function setup_edge_based_rrules(::RedGreenRefinement, topo::UnstructuredGridTop
   if Dc == 2
     faces_list = (ref_nodes,ref_edges,ref_cells)
   else
-    faces_list = (ref_nodes,ref_edges,Int32[],ref_cells)
+    @check length(topo.polytopes) == 1 "Mixed meshes not supported yet."
+    p = first(topo.polytopes)
+    ref_faces = is_simplex(p) ? Int32[] : 1:num_faces(topo,2)
+    faces_list = (ref_nodes,ref_edges,ref_faces,ref_cells)
   end
 
   return rrules, faces_list
@@ -490,18 +497,23 @@ Barycentric refinement method. Each selected cell is refined by adding a new ver
 barycenter of the cell. Then the cell is split into subcells by connecting the new vertex
 to the vertices of the parent cell. Always results in a conformal mesh.
 """
-struct BarycentricRefinementMethod <: EdgeBasedRefinement end
+struct BarycentricRefinement <: EdgeBasedRefinement end
 
 function setup_edge_based_rrules(
-  method::BarycentricRefinementMethod, topo::UnstructuredGridTopology,::Nothing
-)
+  method::BarycentricRefinement, topo::UnstructuredGridTopology{Dc},::Nothing
+) where Dc
   setup_edge_based_rrules(method,topo,collect(1:num_faces(topo,Dc)))
 end
 
 function setup_edge_based_rrules(
-  ::BarycentricRefinementMethod, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}
+  ::BarycentricRefinement, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}
 ) where Dc
-  rrules = lazy_map(BarycentricRefinementRule,CompressedArray(topo.polytopes,topo.cell_type))
+  @check (length(topo.polytopes) == 1) && is_simplex(first(topo.polytopes)) "Only simplex meshes supported"
+
+  ptrs = fill(1,num_faces(topo,Dc))
+  ptrs[cells_to_refine] .= 2
+  p = first(topo.polytopes)
+  rrules = CompressedArray([WhiteRefinementRule(p),BarycentricRefinementRule(p)],ptrs)
 
   ref_nodes = 1:num_faces(topo,0)
   if Dc == 2
@@ -524,7 +536,7 @@ struct RedRefinement <: EdgeBasedRefinementRule end
 struct GreenRefinement{P} <: EdgeBasedRefinementRule end
 struct BlueRefinement{P, Q} <: EdgeBasedRefinementRule end
 struct BlueDoubleRefinement{P} <: EdgeBasedRefinementRule end
-struct BarycentricRefinement <: EdgeBasedRefinementRule end
+struct BarycentricRefinementRule <: EdgeBasedRefinementRule end
 struct PowellSabinRefinement <: EdgeBasedRefinementRule end
 
 _has_interior_point(rr::RefinementRule) = _has_interior_point(rr,RefinementRuleType(rr))
@@ -883,7 +895,7 @@ function BarycentricRefinementRule(p::Polytope)
   reffes = map(x->LagrangianRefFE(Float64,x,1),polys)
 
   ref_grid = UnstructuredDiscreteModel(UnstructuredGrid(coords,conn,reffes,cell_types))
-  return RefinementRule(BarycentricRefinement(),p,ref_grid)
+  return RefinementRule(BarycentricRefinementRule(),p,ref_grid)
 end
 
 function _get_barycentric_refined_faces_list(p::Polytope)
@@ -902,7 +914,7 @@ function _get_barycentric_refined_connectivity(p::Polytope)
     conn_data = [1, 2, 4,
                  2, 3, 4,
                  3, 1, 4]
-    conn_ptrs = [1, 4, 7]
+    conn_ptrs = [1, 4, 7, 10]
     return polys, cell_type, Table(conn_data,conn_ptrs)
   elseif p == TET
     polys     = [TET]
@@ -987,7 +999,7 @@ end
 
 get_relabeled_connectivity(::RefinementRuleType,::RefinementRule,faces_gids) = @notimplemented
 
-function get_relabeled_connectivity(::WithoutRefinement,rr::RefinementRule{<:Polytope{2}},faces_gids)
+function get_relabeled_connectivity(::WithoutRefinement,rr::RefinementRule,faces_gids)
   conn = rr.ref_grid.grid.cell_node_ids
   gids = faces_gids[1]
   new_data = lazy_map(Reindex(gids),conn.data)
@@ -1028,7 +1040,7 @@ function get_relabeled_connectivity(::BlueDoubleRefinement{N},rr::RefinementRule
   return Table(new_data,conn.ptrs)
 end
 
-function get_relabeled_connectivity(::BarycentricRefinement,rr::RefinementRule{<:Polytope{D}},faces_gids) where D
+function get_relabeled_connectivity(::BarycentricRefinementRule,rr::RefinementRule{<:Polytope{D}},faces_gids) where D
   conn = rr.ref_grid.grid.cell_node_ids
   gids = [faces_gids[1]...,faces_gids[D+1]...] # Polytope nodes + Barycenter
   new_data = lazy_map(Reindex(gids),conn.data)
diff --git a/test/AdaptivityTests/EdgeBasedRefinementTests.jl b/test/AdaptivityTests/EdgeBasedRefinementTests.jl
index 9e525318e..b5f8dcce3 100644
--- a/test/AdaptivityTests/EdgeBasedRefinementTests.jl
+++ b/test/AdaptivityTests/EdgeBasedRefinementTests.jl
@@ -104,23 +104,27 @@ end
 
 visualize = false
 
-# Refining meshes of QUADs
+############################################################################################
+### Red-Green refinement
+
+## A) 2D meshes - QUADs 
+
 cart_model = CartesianDiscreteModel((0,1,0,1),(4,4))
 model1     = UnstructuredDiscreteModel(cart_model)
 
-## Homogeneous refinement
+# Homogeneous refinement
 ref_model1 = refine(model1)
 trian1 = Triangulation(ref_model1.model)
 visualize && writevtk(trian1,"test/AdaptivityTests/ref_model1")
 test_grid_transfers(2,model1,ref_model1,1)
 
-## Propagate to all-red
+# Propagate to all-red
 ref_model2 = refine(model1;cells_to_refine=[1,6,11,16])
 trian2 = Triangulation(ref_model2.model)
 visualize && writevtk(trian2,"test/AdaptivityTests/ref_model2")
 test_grid_transfers(2,model1,ref_model2,1)
 
-## Red-Green refinement
+# Red-Green refinement
 ref_model3 = refine(model1;cells_to_refine=[1,6,16])
 trian3 = Triangulation(ref_model3.model)
 visualize && writevtk(trian3,"test/AdaptivityTests/ref_model3")
@@ -131,7 +135,8 @@ trian4 = Triangulation(ref_model4.model)
 visualize && writevtk(trian4,"test/AdaptivityTests/ref_model4")
 #test_grid_transfers(2,model1,ref_model4,1)
 
-# Refining meshes of TRIans
+## B) 2D meshes - TRIs
+
 model2 = simplexify(model1)
 visualize && writevtk(Triangulation(model2),"test/AdaptivityTests/base_model2")
 
@@ -145,21 +150,68 @@ trian6 = Triangulation(ref_model6.model)
 visualize && writevtk(trian6,"test/AdaptivityTests/ref_model6")
 test_grid_transfers(2,model2,ref_model6,1)
 
-# Test newest vertex bisection i.e. longest edge bisection
-# Refine all edges using NVB
-ref_model7 = refine(model2, refinement_method = "nvb")
+## C) 3D meshes - HEXs
+
+cart_model = CartesianDiscreteModel((0,1,0,1,0,1),(2,2,2))
+model3     = UnstructuredDiscreteModel(cart_model)
+
+ref_model7 = refine(model3)
 trian7 = Triangulation(ref_model7.model)
-visualize && writevtk(trian7, "test/AdaptivityTests/ref_model7")
-test_grid_transfers(2, model2, ref_model7, 1)
+visualize && writevtk(trian7,"test/AdaptivityTests/ref_model7")
+test_grid_transfers(3,model3,ref_model7,1)
+
+## D) 3D meshes - TETs
+
+model4 = simplexify(model3)
+ref_model8 = refine(model4)
+trian8 = Triangulation(ref_model8.model)
+visualize && writevtk(trian8,"test/AdaptivityTests/ref_model8")
+test_grid_transfers(3,model4,ref_model8,1)
+
+############################################################################################
+### Newest Vertex Bisection refinement (longest edge bisection)
 
 # Refine all edges using NVB
-# Mark edges such that blue, and double_blue refinement are triggered
-ref_model8 = refine(model2, refinement_method = "nvb", cells_to_refine = [4, 9])
-trian8 = Triangulation(ref_model8)
-visualize && writevtk(trian8, "test/AdaptivityTests/ref_model8")
-ref_model9 = refine(ref_model8, refinement_method = "nvb", cells_to_refine = [1, 3, 4, 11])
+ref_model9 = refine(model2, refinement_method = "nvb")
 trian9 = Triangulation(ref_model9.model)
-visualize && writevtk(trian9, "test/AdaptivityTests/ref_model9")
-test_grid_transfers(2, ref_model8, ref_model9, 1)
+visualize && writevtk(trian9, "test/AdaptivityTests/ref_model7")
+test_grid_transfers(2, model2, ref_model9, 1)
+
+# Refine all edges using NVB
+# Mark edges such that blue, and double_blue refinement are triggered
+ref_model10 = refine(model2, refinement_method = "nvb", cells_to_refine = [4, 9])
+trian10 = Triangulation(ref_model10)
+visualize && writevtk(trian10, "test/AdaptivityTests/ref_model8")
+ref_model11 = refine(ref_model10, refinement_method = "nvb", cells_to_refine = [1, 3, 4, 11])
+trian11 = Triangulation(ref_model11.model)
+visualize && writevtk(trian11, "test/AdaptivityTests/ref_model9")
+test_grid_transfers(2, ref_model10, ref_model11, 1)
+
+############################################################################################
+### Barycentric refinement
+
+## A) 2D meshes - TRIs
+
+ref_model12 = refine(model2, refinement_method = "barycentric")
+trian12 = Triangulation(ref_model12.model)
+visualize && writevtk(trian12, "test/AdaptivityTests/ref_model12")
+test_grid_transfers(2, model2, ref_model12, 1)
+
+ref_model13 = refine(model2, refinement_method = "barycentric", cells_to_refine = [1, 6, 8])
+trian13 = Triangulation(ref_model13.model)
+visualize && writevtk(trian13, "test/AdaptivityTests/ref_model13")
+test_grid_transfers(2, model2, ref_model13, 1)
+
+## B) 3D meshes - TETs
+
+ref_model14 = refine(model4, refinement_method = "barycentric")
+trian14 = Triangulation(ref_model14.model)
+visualize && writevtk(trian13, "test/AdaptivityTests/ref_model13")
+test_grid_transfers(3, model4, ref_model14, 1)
+
+ref_model15 = refine(model4, refinement_method = "barycentric", cells_to_refine = [1, 6, 8])
+trian15 = Triangulation(ref_model15.model)
+visualize && writevtk(trian15, "test/AdaptivityTests/ref_model15")
+test_grid_transfers(3, model4, ref_model15, 1)
 
 end

From 4f6725950d3cda74ba7c593cf02dd728dc5715dd Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Tue, 9 Jul 2024 00:27:27 +1000
Subject: [PATCH 5/8] Bugfix in glues when not all children are in partition

---
 src/Adaptivity/AdaptivityGlues.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/Adaptivity/AdaptivityGlues.jl b/src/Adaptivity/AdaptivityGlues.jl
index bfb434169..ed8fc4ea5 100644
--- a/src/Adaptivity/AdaptivityGlues.jl
+++ b/src/Adaptivity/AdaptivityGlues.jl
@@ -219,6 +219,7 @@ function get_d_to_fface_to_cface(glue::AdaptivityGlue{<:RefinementGlue},
   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)
+  fcell_to_child_id = glue.n2o_cell_to_child_id
 
   # Global data, concerning the complete meshes
   ccell_to_fcell = glue.o2n_faces_map
@@ -234,7 +235,8 @@ function get_d_to_fface_to_cface(glue::AdaptivityGlue{<:RefinementGlue},
       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 fcell in ccell_to_fcell[ccell]
+      child_id = fcell_to_child_id[fcell]
       # 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
@@ -244,7 +246,7 @@ function get_d_to_fface_to_cface(glue::AdaptivityGlue{<:RefinementGlue},
           # 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]
+          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]

From 8c3c55ac2da51289f28a8f533494e2953d1b30eb Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Tue, 9 Jul 2024 17:34:18 +1000
Subject: [PATCH 6/8] Minor

---
 src/Adaptivity/EdgeBasedRefinement.jl | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl
index a60f33b91..fcdf9c761 100644
--- a/src/Adaptivity/EdgeBasedRefinement.jl
+++ b/src/Adaptivity/EdgeBasedRefinement.jl
@@ -195,10 +195,6 @@ function get_refined_cell_to_vertex_map(
     return (N,E,F,C)
   end
 
-  display(faces_list)
-  display(d_to_faces_reindexing)
-  display(d_to_offset)
-
   # Allocate ptr and data arrays for new connectivity
   nC_new    = sum(rr -> num_subcells(rr), rrules)
   nData_new = sum(rr -> sum(length,rr.ref_grid.grid.cell_node_ids), rrules)

From 6950527e069274b443fbb9f5607d03393f70eb07 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 11 Jul 2024 15:07:43 +1000
Subject: [PATCH 7/8] Added Simplexify refinement method

---
 src/Adaptivity/AdaptedDiscreteModels.jl       |  1 +
 src/Adaptivity/Adaptivity.jl                  |  1 +
 src/Adaptivity/AdaptivityGlues.jl             | 28 +++++++++++++++++
 src/Adaptivity/EdgeBasedRefinement.jl         | 26 +---------------
 src/Adaptivity/SimplexifyRefinement.jl        | 30 +++++++++++++++++++
 .../EdgeBasedRefinementTests.jl               | 15 ++++++++++
 6 files changed, 76 insertions(+), 25 deletions(-)
 create mode 100644 src/Adaptivity/SimplexifyRefinement.jl

diff --git a/src/Adaptivity/AdaptedDiscreteModels.jl b/src/Adaptivity/AdaptedDiscreteModels.jl
index 9c23d81b2..d70c0273e 100644
--- a/src/Adaptivity/AdaptedDiscreteModels.jl
+++ b/src/Adaptivity/AdaptedDiscreteModels.jl
@@ -102,6 +102,7 @@ function string_to_refinement(refinement_method::String, model)
   refinement_method == "red_green" && return RedGreenRefinement()
   refinement_method == "nvb" && return NVBRefinement(model)
   refinement_method == "barycentric" && return BarycentricRefinement()
+  refinement_method == "simplexify" && return SimplexifyRefinement()
   error("refinement_method $refinement_method not recognized")
 end
 
diff --git a/src/Adaptivity/Adaptivity.jl b/src/Adaptivity/Adaptivity.jl
index f16d80d8e..2a9b180c1 100644
--- a/src/Adaptivity/Adaptivity.jl
+++ b/src/Adaptivity/Adaptivity.jl
@@ -49,5 +49,6 @@ include("AdaptedDiscreteModels.jl")
 include("AdaptedTriangulations.jl")
 include("CompositeQuadratures.jl")
 include("EdgeBasedRefinement.jl")
+include("SimplexifyRefinement.jl")
 
 end # module
diff --git a/src/Adaptivity/AdaptivityGlues.jl b/src/Adaptivity/AdaptivityGlues.jl
index ed8fc4ea5..3c1cfc2c1 100644
--- a/src/Adaptivity/AdaptivityGlues.jl
+++ b/src/Adaptivity/AdaptivityGlues.jl
@@ -297,4 +297,32 @@ function refine_face_labeling(coarse_labeling::FaceLabeling,
   end
   
   return Geometry.FaceLabeling(d_to_dface_to_entity,tag_to_entities,tag_to_name)  
+end
+
+"""
+    blocked_refinement_glue(rrules::AbstractVector{<:RefinementRule})
+
+Given an array of RefinementRules for each coarse cell, returns an AdaptivityGlue
+where children from the same parent are placed in contiguous blocks.
+"""
+function blocked_refinement_glue(
+  rrules::AbstractVector{<:RefinementRule{<:Polytope{Dc}}}
+) where Dc
+  nC_old = length(rrules)
+  nC_new = sum(num_subcells,rrules)
+
+  f2c_cell_map      = Vector{Int}(undef,nC_new)
+  fcell_to_child_id = Vector{Int}(undef,nC_new)
+
+  k = 1
+  for iC = 1:nC_old
+    rr = rrules[iC]
+    range = k:k+num_subcells(rr)-1
+    f2c_cell_map[range] .= iC
+    fcell_to_child_id[range] .= collect(1:num_subcells(rr))
+    k += num_subcells(rr)
+  end
+
+  f2c_faces_map = [(d==Dc) ? f2c_cell_map : Int[] for d in 0:Dc]
+  return AdaptivityGlue(f2c_faces_map,fcell_to_child_id,rrules)
 end
\ No newline at end of file
diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl
index fcdf9c761..6039c4dd7 100644
--- a/src/Adaptivity/EdgeBasedRefinement.jl
+++ b/src/Adaptivity/EdgeBasedRefinement.jl
@@ -63,7 +63,7 @@ function refine(method::EdgeBasedRefinement,model::UnstructuredDiscreteModel{Dc,
     get_cell_type(topo),
     OrientationStyle(topo)
   )
-  glue = get_edge_based_refinement_glue(topo,model.grid_topology,rrules)
+  glue = blocked_refinement_glue(rrules)
   labels = refine_face_labeling(coarse_labels,glue,model.grid_topology,topo)
   ref_model = UnstructuredDiscreteModel(grid,topo,labels)
   return AdaptedDiscreteModel(ref_model,model,glue)
@@ -82,30 +82,6 @@ function refine_edge_based_topology(
   return UnstructuredGridTopology(coords_new,c2n_map_new,cell_type_new,polys_new,orientation)
 end
 
-function get_edge_based_refinement_glue(
-  ftopo ::UnstructuredGridTopology{Dc},
-  ctopo ::UnstructuredGridTopology{Dc},
-  rrules::AbstractVector{<:RefinementRule}
-) where {Dc}
-  nC_old = num_faces(ctopo,Dc)
-  nC_new = num_faces(ftopo,Dc)
-
-  f2c_cell_map      = Vector{Int}(undef,nC_new)
-  fcell_to_child_id = Vector{Int}(undef,nC_new)
-
-  k = 1
-  for iC = 1:nC_old
-    rr = rrules[iC]
-    range = k:k+num_subcells(rr)-1
-    f2c_cell_map[range] .= iC
-    fcell_to_child_id[range] .= collect(1:num_subcells(rr))
-    k += num_subcells(rr)
-  end
-
-  f2c_faces_map = [(d==Dc) ? f2c_cell_map : Int[] for d in 0:Dc]
-  return AdaptivityGlue(f2c_faces_map,fcell_to_child_id,rrules)
-end
-
 """
     setup_edge_based_rrules(method::EdgeBasedRefinement,topo::UnstructuredGridTopology,cells_to_refine)
 
diff --git a/src/Adaptivity/SimplexifyRefinement.jl b/src/Adaptivity/SimplexifyRefinement.jl
new file mode 100644
index 000000000..1239b5e67
--- /dev/null
+++ b/src/Adaptivity/SimplexifyRefinement.jl
@@ -0,0 +1,30 @@
+
+"""
+    struct SimplexifyRefinement <: AdaptivityMethod
+
+Equivalent to `simplexify`, but keeps track of the parent-child relationship between 
+the original and the refined model.
+"""
+struct SimplexifyRefinement <: AdaptivityMethod end
+
+function refine(::SimplexifyRefinement,model::UnstructuredDiscreteModel{Dc,Dp}) where {Dc,Dp}
+  ref_model = simplexify(model)
+
+  polys = get_polytopes(model)
+  ctype = get_cell_type(model)
+  rrules = expand_cell_data(map(p -> SimplexifyRefinementRule(p),polys),ctype)
+  glue = blocked_refinement_glue(rrules)
+  
+  return AdaptedDiscreteModel(ref_model,model,glue)
+end
+
+struct SimplexifyRefinementRule <: RefinementRuleType end
+
+function SimplexifyRefinementRule(poly::Polytope)
+  conn, sp = simplexify(poly)
+  coords = get_vertex_coordinates(poly)
+  reffes = [LagrangianRefFE(Float64,sp,1)]
+  cell_types = fill(1,length(conn))
+  ref_grid = UnstructuredDiscreteModel(UnstructuredGrid(coords,Table(conn),reffes,cell_types))
+  return RefinementRule(SimplexifyRefinementRule(),poly,ref_grid)
+end
diff --git a/test/AdaptivityTests/EdgeBasedRefinementTests.jl b/test/AdaptivityTests/EdgeBasedRefinementTests.jl
index b5f8dcce3..4ed3f9c2f 100644
--- a/test/AdaptivityTests/EdgeBasedRefinementTests.jl
+++ b/test/AdaptivityTests/EdgeBasedRefinementTests.jl
@@ -214,4 +214,19 @@ trian15 = Triangulation(ref_model15.model)
 visualize && writevtk(trian15, "test/AdaptivityTests/ref_model15")
 test_grid_transfers(3, model4, ref_model15, 1)
 
+############################################################################################
+### Simplexify refinement
+
+## A) 2D meshes - QUADs
+ref_model16 = refine(model1, refinement_method = "simplexify")
+trian16 = Triangulation(ref_model16.model)
+visualize && writevtk(trian16, "test/AdaptivityTests/ref_model16")
+test_grid_transfers(2, model1, ref_model16, 1)
+
+## B) 3D meshes - HEXs
+ref_model17 = refine(model3, refinement_method = "simplexify")
+trian17 = Triangulation(ref_model17.model)
+visualize && writevtk(trian17, "test/AdaptivityTests/ref_model17")
+test_grid_transfers(3, model3, ref_model17, 1)
+
 end

From c8e13a5e29f15eaca60708cd5c9010c666d40e83 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 11 Jul 2024 17:08:59 +1000
Subject: [PATCH 8/8] Updated NEWS and bumped release

---
 NEWS.md      | 6 ++++--
 Project.toml | 2 +-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 0298fb970..5de323b2d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,13 +5,15 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased]
+## [0.18.3] - 2024-07-11
 
 ### Added
 
+- Added more features to Adaptivity. Notably: 3D uniform edge-based refinement for HEX and TETs. Barycentric refinement for simplices. Simplexify as a new refinement strategy. Since PR[#1013](https://github.com/gridap/Gridap.jl/pull/1013).
+
 - Define `GeneralPolytope` that represents general polytopes in 2 and 3 dimensions. Since PR[#1006](https://github.com/gridap/Gridap.jl/pull/1006).
 
-## [0.18.2] - 2024-05-02 
+## [0.18.2] - 2024-05-02
 
 ### Fixed
 
diff --git a/Project.toml b/Project.toml
index 3ee248317..1aae37493 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 authors = ["Santiago Badia <santiago.badia@monash.edu>", "Francesc Verdugo <f.verdugo.rojano@vu.nl>", "Alberto F. Martin <alberto.f.martin@anu.edu.au>"]
-version = "0.18.2"
+version = "0.18.3"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"