diff --git a/src/Adaptivity.jl b/src/Adaptivity.jl
index c3085ba..1d9ef70 100644
--- a/src/Adaptivity.jl
+++ b/src/Adaptivity.jl
@@ -35,6 +35,25 @@ function Adaptivity.get_adaptivity_glue(model::DistributedAdaptedDiscreteModel)
   return map(Adaptivity.get_adaptivity_glue,local_views(model))
 end
 
+function Adaptivity.refine(
+  cmodel::DistributedAdaptedDiscreteModel{Dc},args...;kwargs...
+) where Dc
+  # Local cmodels are AdaptedDiscreteModels. To correctly dispatch, we need to
+  # extract the underlying models, then refine.
+  _cmodel = get_model(cmodel)
+  _fmodel = refine(_cmodel,args...;kwargs...)
+
+  # Now the issue is that the local parents are not pointing to local_views(cmodel).
+  # We have to fix that...
+  fmodel = GenericDistributedDiscreteModel(
+    map(get_model,local_views(_fmodel)),
+    get_cell_gids(_fmodel);
+    metadata=_fmodel.metadata
+  )
+  glues = get_adaptivity_glue(_fmodel)
+  return DistributedAdaptedDiscreteModel(fmodel,cmodel,glues)
+end
+
 # RedistributeGlue : Redistributing discrete models
 
 """
@@ -494,27 +513,6 @@ function Adaptivity.refine(
   Adaptivity.refine(cmodel,Tuple(fill(refs,Dc)))
 end
 
-function Adaptivity.refine(
-  cmodel::DistributedAdaptedDiscreteModel{Dc},
-  refs::NTuple{Dc,<:Integer}
-) where Dc
-
-  # Local cmodels are AdaptedDiscreteModels. To correctly dispatch, we need to
-  # extract the underlying models, then refine.
-  _cmodel = get_model(cmodel)
-  _fmodel = refine(_cmodel,refs)
-
-  # Now the issue is that the local parents are not pointing to local_views(cmodel).
-  # We have to fix that...
-  fmodel = GenericDistributedDiscreteModel(
-    map(get_model,local_views(_fmodel)),
-    get_cell_gids(_fmodel);
-    metadata=_fmodel.metadata
-  )
-  glues = get_adaptivity_glue(_fmodel)
-  return DistributedAdaptedDiscreteModel(fmodel,cmodel,glues)
-end
-
 function Adaptivity.refine(
   cmodel::DistributedCartesianDiscreteModel{Dc},
   refs::NTuple{Dc,<:Integer},
@@ -718,8 +716,8 @@ function get_cartesian_redistribute_glue(
       #   - I don't send anything. 
       #   - I receive all my owned cells in the new model.
       old2new = Int[]
-      new2old = fill(0,num_cells(new_model))
-      ids_rcv = collect(own_to_local(new_ids))
+      new2old = fill(zero(Int),num_cells(new_model))
+      ids_rcv = collect(Int,own_to_local(new_ids))
       ids_snd = Int[]
     end
 
@@ -793,7 +791,7 @@ function get_redistributed_face_labeling(
         old_face2entity = old_labels.d_to_dface_to_entity[Df+1]
 
         old_cell_to_face_entity = Table(
-          lazy_map(Reindex(old_face2entity),old_cell2face.data),
+          collect(Int32,lazy_map(Reindex(old_face2entity),old_cell2face.data)), # Avoid if possible
           old_cell2face.ptrs
         )
       else
@@ -868,3 +866,250 @@ function get_redistributed_face_labeling(
   new_labels = map(FaceLabeling,new_d_to_dface_to_entity,new_tag_to_entities,new_tag_to_name)
   return DistributedFaceLabeling(new_labels)
 end
+
+# UnstructuredDiscreteModel refinement
+
+function Adaptivity.refine(
+  cmodel::DistributedUnstructuredDiscreteModel{Dc},args...;kwargs...
+) where Dc
+  fmodels, f_own_to_local = refine_local_models(cmodel,args...;kwargs...)
+  fgids = refine_cell_gids(cmodel,fmodels,f_own_to_local)
+  return GenericDistributedDiscreteModel(fmodels,fgids)
+end
+
+"""
+    refine_coarse_models(cmodel::DistributedDiscreteModel{Dc},args...;kwargs...) where Dc
+
+Given a coarse model, returns the locally refined models. This is done by 
+  - refining the local models serially
+  - filtering out the extra fine layers of ghosts
+We also return the ids of the owned fine cells.
+
+To find the fine cells we want to keep, we have the following criteria: 
+  - Given a fine cell, it is owned iff 
+    A) It's parent cell is owned
+  - Given a fine cell, it is a ghost iff not(A) and 
+    B) It has at least one neighbor with a non-ghost parent
+
+Instead of checking A and B, we do the following: 
+  - We mark fine owned cells by checking A 
+  - If a cell is owned, we set it's fine neighbors as owned or ghost
+"""
+function refine_local_models(
+  cmodel::DistributedDiscreteModel{Dc},args...;kwargs...
+) where Dc
+  cgids = partition(get_cell_gids(cmodel))
+  cmodels = local_views(cmodel)
+
+  # Refine models locally
+  fmodels = map(cmodels) do cmodel
+    refine(cmodel,args...;kwargs...)
+  end
+
+  # Select fine cells we want to keep
+  Df = Dc-1 # Dimension used to find neighboring cells\
+  f_own_or_ghost_ids, f_own_ids = map(cgids,cmodels,fmodels) do cgids,cmodel,fmodel
+    glue = get_adaptivity_glue(fmodel)
+    f2c_map = glue.n2o_faces_map[Dc+1]
+    child_map = glue.n2o_cell_to_child_id
+
+    ftopo = get_grid_topology(fmodel)
+    f_cell_to_facet = Geometry.get_faces(ftopo,Dc,Df)
+    f_facet_to_cell = Geometry.get_faces(ftopo,Df,Dc)
+    f_cell_to_facet_cache = array_cache(f_cell_to_facet)
+    f_facet_to_cell_cache = array_cache(f_facet_to_cell)
+    c_l2o_map = local_to_own(cgids)
+    
+    f_own_mask = fill(false,length(f2c_map))
+    f_own_or_ghost_mask = fill(false,length(f2c_map))
+    for (fcell,ccell) in enumerate(f2c_map)
+      if !iszero(c_l2o_map[ccell])
+        f_own_mask[fcell] = true
+        facets = getindex!(f_cell_to_facet_cache,f_cell_to_facet,fcell)
+        for facet in facets
+          facet_cells = getindex!(f_facet_to_cell_cache,f_facet_to_cell,facet)
+          for facet_cell in facet_cells
+            f_own_or_ghost_mask[facet_cell] = true
+          end
+        end
+      end
+    end
+
+    f_own_or_ghost_ids = findall(f_own_or_ghost_mask)
+    f_own_ids = findall(i -> f_own_mask[i],f_own_or_ghost_ids) # ModelPortion numeration
+
+    return f_own_or_ghost_ids, f_own_ids
+  end |> tuple_of_arrays
+
+  # Filter out local models
+  filtered_fmodels = map(fmodels,f_own_or_ghost_ids) do fmodel,f_own_or_ghost_ids
+    model = DiscreteModelPortion(get_model(fmodel),f_own_or_ghost_ids).model
+    parent = get_parent(fmodel)
+
+    _glue = get_adaptivity_glue(fmodel)
+    n2o_faces_map = Vector{Vector{Int}}(undef,Dc+1)
+    n2o_faces_map[Dc+1] = _glue.n2o_faces_map[Dc+1][f_own_or_ghost_ids]
+    n2o_cell_to_child_id = _glue.n2o_cell_to_child_id[f_own_or_ghost_ids]
+    rrules = _glue.refinement_rules
+    glue = AdaptivityGlue(n2o_faces_map,n2o_cell_to_child_id,rrules)
+    return AdaptedDiscreteModel(model,parent,glue)
+  end
+
+  return filtered_fmodels, f_own_ids
+end
+
+"""
+    refine_cell_gids(
+      cmodel::DistributedDiscreteModel{Dc},
+      fmodels::AbstractArray{<:DiscreteModel{Dc}}
+    ) where Dc
+
+Given a coarse model and it's local refined models, returns the gids of the fine model.
+The gids are computed as follows: 
+  - First, we create a global numbering for the owned cells by adding an owner-based offset to the local 
+    cell ids (such that cells belonging to the first processor are numbered first). This is 
+    quite standard.
+  - The complicated part is making this numeration consistent, i.e communicating gids of the 
+    ghost cells. To do so, each processor selects it's ghost fine cells, and requests their 
+    global ids by sending two keys:
+      1. The global id of the coarse parent
+      2. The child id of the fine cell
+"""
+function refine_cell_gids(
+  cmodel::DistributedDiscreteModel{Dc},
+  fmodels::AbstractArray{<:DiscreteModel{Dc}},
+  f_own_to_local::AbstractArray{<:AbstractArray{Int}},
+) where Dc
+
+  cgids = partition(get_cell_gids(cmodel))
+  cmodels = local_views(cmodel)
+  ranks = linear_indices(cgids)
+
+  # Create own numbering (without ghosts)
+  num_f_owned_cells = map(length,f_own_to_local)
+  num_f_gids = reduce(+,num_f_owned_cells)
+  first_f_gid = scan(+,num_f_owned_cells,type=:exclusive,init=1)
+  
+  own_fgids = map(ranks,first_f_gid,num_f_owned_cells) do rank,first_f_gid,num_f_owned_cells
+    f_o2g = collect(first_f_gid:first_f_gid+num_f_owned_cells-1)
+    own   = OwnIndices(num_f_gids,rank,f_o2g)
+    ghost = GhostIndices(num_f_gids) # No ghosts
+    return OwnAndGhostIndices(own,ghost)
+  end
+  
+  # Select ghost fine cells local ids
+  parts_rcv, parts_snd = assembly_neighbors(cgids);
+  lids_snd = map(parts_snd,cgids,fmodels) do parts_snd,cgids,fmodel
+    glue = get_adaptivity_glue(fmodel)
+    f2c_map = glue.n2o_faces_map[Dc+1]
+    c_owners = local_to_owner(cgids)
+    return JaggedArray(map(p -> findall(parent -> c_owners[parent] == p,f2c_map),parts_snd))
+  end
+  
+  # Given the local ids of the fine cells we want to get info on, we 
+  # collect two keys: 
+  #   1. The global id of the coarse parent
+  #   2. The child id of the fine cell
+  parent_gids_snd, child_ids_snd = map(cgids,cmodels,fmodels,lids_snd) do cgids,cmodel,fmodel,lids_snd
+    glue = get_adaptivity_glue(fmodel)
+    f2c_map = glue.n2o_faces_map[Dc+1]
+    child_map = glue.n2o_cell_to_child_id
+  
+    c_l2g_map = local_to_global(cgids)
+    c_owners  = local_to_owner(cgids)
+  
+    parent_gids, child_ids, owners = map(lids_snd.data) do fcell
+      ccell = f2c_map[fcell]
+      return c_l2g_map[ccell], child_map[fcell], c_owners[ccell]
+    end |> tuple_of_arrays
+  
+    ptrs = lids_snd.ptrs
+    return JaggedArray(parent_gids,ptrs), JaggedArray(child_ids,ptrs)
+  end |> tuple_of_arrays;
+  
+  # We exchange the keys
+  parts_rcv, parts_snd = assembly_neighbors(cgids);
+  graph = ExchangeGraph(parts_snd,parts_rcv)
+  t1 = exchange(parent_gids_snd,graph)
+  t2 = exchange(child_ids_snd,graph)
+  parent_gids_rcv = fetch(t1)
+  child_ids_rcv = fetch(t2)
+  
+  # We process the received keys, and collect the global ids of the fine cells
+  # that have been requested by our neighbors.
+  child_gids_rcv = map(
+    cgids,own_fgids,f_own_to_local,cmodels,fmodels,parent_gids_rcv,child_ids_rcv
+  ) do cgids,own_fgids,f_own_to_local,cmodel,fmodel,parent_gids_rcv,child_ids_rcv
+    glue = get_adaptivity_glue(fmodel)
+    c2f_map = glue.o2n_faces_map
+    child_map = glue.n2o_cell_to_child_id
+    c2f_map_cache = array_cache(c2f_map)
+  
+    parent_lids = to_local!(parent_gids_rcv.data,cgids)
+    child_ids = child_ids_rcv.data
+  
+    f_local_to_own = Arrays.find_inverse_index_map(f_own_to_local)
+  
+    child_lids = map(parent_lids,child_ids) do ccell, child_id
+      fcells = getindex!(c2f_map_cache,c2f_map,ccell)
+      pos = findfirst(fcell -> child_map[fcell] == child_id, fcells)
+      return f_local_to_own[fcells[pos]]
+    end
+  
+    child_gids = to_global!(child_lids,own_fgids)
+    return JaggedArray(child_gids,parent_gids_rcv.ptrs)
+  end
+  
+  # We exchange back the information
+  graph = ExchangeGraph(parts_rcv,parts_snd)
+  t = exchange(child_gids_rcv,graph)
+  child_gids_snd = fetch(t)
+  
+  # We finally can create the global numeration of the fine cells by piecing together: 
+  #   1. The (local ids,global ids) of the owned fine cells
+  #   2. The (owners,local ids,global ids) of the ghost fine cells
+  fgids = map(ranks,f_own_to_local,own_fgids,lids_snd,child_gids_snd) do rank, own_lids, own_gids, ghost_lids, ghost_gids
+    own2global = own_to_global(own_gids)
+  
+    n_nbors = length(ghost_lids)
+    n_own   = length(own_lids)
+    n_ghost = length(ghost_lids.data)
+    local2global = fill(0,n_own+n_ghost)
+    local2owner = fill(0,n_own+n_ghost)
+  
+    # Own cells
+    for (oid,lid) in enumerate(own_lids)
+      local2global[lid] = own2global[oid]
+      local2owner[lid]  = rank
+    end
+    
+    # Ghost cells
+    for nbor in 1:n_nbors
+      for i in ghost_lids.ptrs[nbor]:ghost_lids.ptrs[nbor+1]-1
+        lid = ghost_lids.data[i]
+        gid = ghost_gids.data[i]
+        local2global[lid] = gid
+        local2owner[lid]  = nbor
+      end
+    end
+    return LocalIndices(num_f_gids,rank,local2global,local2owner)
+  end
+
+  return PRange(fgids)
+end
+
+function refine_cell_gids(
+  cmodel::DistributedDiscreteModel{Dc},
+  fmodels::AbstractArray{<:DiscreteModel{Dc}}
+) where Dc
+  cgids = get_cell_gids(cmodel)
+  f_own_to_local = map(cgids,fmodels) do cgids,fmodel
+    glue = get_adaptivity_glue(fmodel)
+    f2c_map = glue.n2o_faces_map[Dc+1]
+    @assert isa(f2c_map,Vector) "Only uniform refinement is supported!"
+  
+    c_l2o_map = local_to_own(cgids)
+    return findall(parent -> !iszero(c_l2o_map[parent]),f2c_map)
+  end
+  return refine_cell_gids(cmodel,fmodels,f_own_to_local)
+end
diff --git a/src/Geometry.jl b/src/Geometry.jl
index 3263e42..57686eb 100644
--- a/src/Geometry.jl
+++ b/src/Geometry.jl
@@ -415,6 +415,18 @@ function Geometry.DiscreteModel(
   GenericDistributedDiscreteModel(models,gids)
 end
 
+# UnstructuredDiscreteModel
+
+const DistributedUnstructuredDiscreteModel{Dc,Dp,A,B,C} = 
+  GenericDistributedDiscreteModel{Dc,Dp,<:AbstractArray{<:UnstructuredDiscreteModel},B,C}
+
+function Geometry.UnstructuredDiscreteModel(model::GenericDistributedDiscreteModel)
+  return GenericDistributedDiscreteModel(
+    map(UnstructuredDiscreteModel,local_views(model)),
+    get_cell_gids(model),
+  )
+end
+
 # Triangulation
 
 # We do not inherit from Triangulation on purpose.
diff --git a/test/AdaptivityUnstructuredTests.jl b/test/AdaptivityUnstructuredTests.jl
index 0e19f58..22c415f 100644
--- a/test/AdaptivityUnstructuredTests.jl
+++ b/test/AdaptivityUnstructuredTests.jl
@@ -1,4 +1,4 @@
-
+using FillArrays
 using PartitionedArrays, GridapDistributed
 using Gridap
 using Gridap.Adaptivity, Gridap.Arrays, Gridap.Geometry
@@ -9,171 +9,10 @@ ranks = with_debug() do distribute
   distribute(LinearIndices((2,)))
 end
 
-Dc = 2
-
 _cmodel = CartesianDiscreteModel(ranks,(2,1),(0,1,0,1),(4,4))
-cmodel = GenericDistributedDiscreteModel(
-  map(UnstructuredDiscreteModel,local_views(_cmodel)),
-  get_cell_gids(_cmodel),
-)
-
-# A) Refine local models
-cgids = partition(get_cell_gids(cmodel))
-cmodels = local_views(cmodel)
-fmodels = map(cmodels) do cmodel
-  refine(cmodel)
-end
-
-# B) Remove extra ghost cells
-
-f_ghost_lids = map(cgids,cmodels,fmodels) do cgids,cmodel,fmodel
-  glue = get_adaptivity_glue(fmodel)
-  f2c_map = glue.n2o_faces_map[Dc+1]
-  child_map = glue.n2o_cell_to_child_id
-
-  ftopo = get_grid_topology(fmodel)
-  f_cell_to_cell = Geometry.get_faces(ftopo,Dc,Dc) # Not what we want, we need to go through faces
-  f_cell_to_cell_cache = array_cache(f_cell_to_cell)
-  c_l2o_map = local_to_own(cgids)
-  
-  f_ghost_mask = fill(false,length(f2c_map))
-  for (fcell,ccell) in enumerate(f2c_map)
-    fine_nbors = getindex!(f_cell_to_cell_cache,f_cell_to_cell,fcell)
-    println(fine_nbors)
-    A = iszero(c_l2o_map[ccell]) # Parent is ghost
-    B = any(nbor -> !iszero(c_l2o_map[f2c_map[nbor]]), fine_nbors) # Has neighbor with non-ghost parent
-    f_ghost_mask[fcell] = A && B
-  end
-
-  return findall(f_ghost_mask)
-end
-
-
-# C) Create global numeration of own cells by
-#   1. Counting the number of owned fine cells in each model (children of owned coarse cells)
-#   2. Computing the first global id of each model, and the maximum gobal id
-#   3. Creating a global numeration of owned fine cells
-f_own_to_local = map(cgids,cmodels,fmodels) do cgids,cmodel,fmodel
-  glue = get_adaptivity_glue(fmodel)
-  f2c_map = glue.n2o_faces_map[Dc+1]
-  child_map = glue.n2o_cell_to_child_id
-  @assert isa(f2c_map,Vector) "Only uniform refinement is supported!"
-
-  c_l2o_map = local_to_own(cgids)
-  return findall(parent -> !iszero(c_l2o_map[parent]),f2c_map)
-end
-
-num_f_owned_cells = map(length,f_own_to_local)
-num_f_gids = reduce(+,num_f_owned_cells)
-first_f_gid = scan(+,num_f_owned_cells,type=:exclusive,init=1)
-
-own_fgids = map(ranks,first_f_gid,num_f_owned_cells) do rank,first_f_gid,num_f_owned_cells
-  f_o2g = collect(first_f_gid:first_f_gid+num_f_owned_cells-1)
-  own   = OwnIndices(num_f_gids,rank,f_o2g)
-  ghost = GhostIndices(num_f_gids) # No ghosts
-  return OwnAndGhostIndices(own,ghost)
-end
+cmodel = UnstructuredDiscreteModel(_cmodel)
 
-# D) Exchange ghost gids by sending two keys: 
-#   1. The global id of the coarse parent
-#   2. The child id of the fine cell
+fmodel = refine(cmodel)
 
-# This has extra ghosts, but will do for now (see above)
-parts_rcv, parts_snd = assembly_neighbors(cgids);
-lids_snd = map(ranks,parts_snd,cgids,cmodels,fmodels) do r,parts_snd,cgids,cmodel,fmodel
-  glue = get_adaptivity_glue(fmodel)
-  f2c_map = glue.n2o_faces_map[Dc+1]
-  child_map = glue.n2o_cell_to_child_id
-  c_owners = local_to_owner(cgids)
-  return JaggedArray(map(p -> findall(parent -> c_owners[parent] == p,f2c_map),parts_snd))
-end
-
-# Given the local ids of the fine cells we want to get info on, we 
-# collect two keys: 
-#   1. The global id of the coarse parent
-#   2. The child id of the fine cell
-parent_gids_snd, child_ids_snd = map(cgids,cmodels,fmodels,lids_snd) do cgids,cmodel,fmodel,lids_snd
-  glue = get_adaptivity_glue(fmodel)
-  f2c_map = glue.n2o_faces_map[Dc+1]
-  child_map = glue.n2o_cell_to_child_id
-
-  c_l2g_map = local_to_global(cgids)
-  c_owners  = local_to_owner(cgids)
-
-  parent_gids, child_ids, owners = map(lids_snd.data) do fcell
-    ccell = f2c_map[fcell]
-    return c_l2g_map[ccell], child_map[fcell], c_owners[ccell]
-  end |> tuple_of_arrays
-
-  ptrs = lids_snd.ptrs
-  return JaggedArray(parent_gids,ptrs), JaggedArray(child_ids,ptrs)
-end |> tuple_of_arrays;
-
-# We exchange the keys
-parts_rcv, parts_snd = assembly_neighbors(cgids);
-graph = ExchangeGraph(parts_snd,parts_rcv)
-t1 = exchange(parent_gids_snd,graph)
-t2 = exchange(child_ids_snd,graph)
-parent_gids_rcv = fetch(t1)
-child_ids_rcv = fetch(t2)
-
-# We process the received keys, and collect the global ids of the fine cells
-# that have been requested by our neighbors.
-child_gids_rcv = map(
-  cgids,own_fgids,f_own_to_local,cmodels,fmodels,parent_gids_rcv,child_ids_rcv
-) do cgids,own_fgids,f_own_to_local,cmodel,fmodel,parent_gids_rcv,child_ids_rcv
-  glue = get_adaptivity_glue(fmodel)
-  c2f_map = glue.o2n_faces_map
-  child_map = glue.n2o_cell_to_child_id
-  c2f_map_cache = array_cache(c2f_map)
-
-  parent_lids = to_local!(parent_gids_rcv.data,cgids)
-  child_ids = child_ids_rcv.data
 
-  f_local_to_own = Arrays.find_inverse_index_map(f_own_to_local)
-
-  child_lids = map(parent_lids,child_ids) do ccell, child_id
-    fcells = getindex!(c2f_map_cache,c2f_map,ccell)
-    pos = findfirst(fcell -> child_map[fcell] == child_id, fcells)
-    return f_local_to_own[fcells[pos]]
-  end
-
-  child_gids = to_global!(child_lids,own_fgids)
-  return JaggedArray(child_gids,parent_gids_rcv.ptrs)
-end
-
-# We exchange back the information
-graph = ExchangeGraph(parts_rcv,parts_snd)
-t = exchange(child_gids_rcv,graph)
-child_gids_snd = fetch(t)
-
-# We finally can create the global numeration of the fine cells by piecing together: 
-#   1. The (local ids,global ids) of the owned fine cells
-#   2. The (owners,local ids,global ids) of the ghost fine cells
-fgids = map(ranks,f_own_to_local,own_fgids,lids_snd,child_gids_snd) do rank, own_lids, own_gids, ghost_lids, ghost_gids
-  own2global = own_to_global(own_gids)
-
-  n_nbors = length(ghost_lids)
-  n_own   = length(own_lids)
-  n_ghost = length(ghost_lids.data)
-  local2global = fill(0,n_own+n_ghost)
-  local2owner = fill(0,n_own+n_ghost)
-
-  # Own cells
-  for (oid,lid) in enumerate(own_lids)
-    local2global[lid] = own2global[oid]
-    local2owner[lid]  = rank
-  end
-  
-  # Ghost cells
-  for nbor in 1:n_nbors
-    for i in ghost_lids.ptrs[nbor]:ghost_lids.ptrs[nbor+1]-1
-      lid = ghost_lids.data[i]
-      gid = ghost_gids.data[i]
-      local2global[lid] = gid
-      local2owner[lid]  = nbor
-    end
-  end
-  return LocalIndices(num_f_gids,rank,local2global,local2owner)
-end