diff --git a/NEWS.md b/NEWS.md
index 7d7bb3a..325955f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 - Added support for distributed level-set geometries. Since PR[#99](https://github.com/gridap/GridapEmbedded.jl/pull/99).
+- Refactored the distributed code to allow for ghosted/unghosted geometries and triangulations. Since PR[#100](https://github.com/gridap/GridapEmbedded.jl/pull/100).
 
 ## [0.9.5] - 2024-10-18
 
diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl
index ee2ea62..b4cf0b8 100644
--- a/src/Distributed/Distributed.jl
+++ b/src/Distributed/Distributed.jl
@@ -26,19 +26,16 @@ using GridapEmbedded.Interfaces: AbstractEmbeddedDiscretization
 using GridapEmbedded.AgFEM: _touch_aggregated_cells!
 using GridapEmbedded.AgFEM: AggregateCutCellsByThreshold
 using GridapEmbedded.MomentFittedQuadratures: MomentFitted
-using Gridap.Geometry: AppendedTriangulation
+using Gridap.Geometry: AppendedTriangulation, TriangulationView
 using Gridap.Geometry: get_face_to_parent_face
 using Gridap.Arrays: find_inverse_index_map, testitem, return_type
 using Gridap.FESpaces: FESpaceWithLinearConstraints
 using Gridap.FESpaces: _dof_to_DOF, _DOF_to_dof
-using GridapDistributed: DistributedDiscreteModel
-using GridapDistributed: DistributedTriangulation
-using GridapDistributed: DistributedFESpace
-using GridapDistributed: DistributedSingleFieldFESpace
-using GridapDistributed: DistributedMeasure
-using GridapDistributed: add_ghost_cells
-using GridapDistributed: generate_gids
-using GridapDistributed: generate_cell_gids
+
+using GridapDistributed: DistributedDiscreteModel, DistributedTriangulation, DistributedMeasure
+using GridapDistributed: DistributedFESpace, DistributedSingleFieldFESpace
+using GridapDistributed: NoGhost, WithGhost, filter_cells_when_needed, add_ghost_cells
+using GridapDistributed: generate_gids, generate_cell_gids
 using GridapDistributed: _find_vector_type
 
 import GridapEmbedded.AgFEM: aggregate
@@ -48,6 +45,7 @@ import GridapEmbedded.Interfaces: cut_facets
 import GridapEmbedded.Interfaces: EmbeddedBoundary
 import GridapEmbedded.Interfaces: compute_bgfacet_to_inoutcut
 import GridapEmbedded.Interfaces: compute_bgcell_to_inoutcut
+import GridapEmbedded.Interfaces: GhostSkeleton
 import GridapEmbedded.CSG: get_geometry
 import GridapEmbedded.LevelSetCutters: discretize, DiscreteGeometry
 import Gridap.Geometry: Triangulation
@@ -61,6 +59,8 @@ include("DistributedDiscretizations.jl")
 
 include("DistributedDiscreteGeometries.jl")
 
+include("DistributedSubFacetTriangulations.jl")
+
 include("DistributedAgFEM.jl")
 
 include("DistributedQuadratures.jl")
diff --git a/src/Distributed/DistributedAgFEM.jl b/src/Distributed/DistributedAgFEM.jl
index 4501ba7..cdb65ab 100644
--- a/src/Distributed/DistributedAgFEM.jl
+++ b/src/Distributed/DistributedAgFEM.jl
@@ -3,8 +3,8 @@ function AgFEMSpace(
   bgmodel::DistributedDiscreteModel,
   f::DistributedFESpace,
   bgcell_to_bgcellin::AbstractArray{<:AbstractVector},
-  g::DistributedFESpace=f)
-
+  g::DistributedFESpace=f
+)
   bgmodel_gids = get_cell_gids(bgmodel)
   spaces = map(
     local_views(f),
@@ -13,9 +13,7 @@ function AgFEMSpace(
     local_views(bgmodel_gids)) do f,bgcell_to_bgcellin,g,gids
       AgFEMSpace(f,bgcell_to_bgcellin,g,local_to_global(gids))
   end
-  trians = map(get_triangulation,local_views(f))
-  trian = DistributedTriangulation(trians,bgmodel)
-  trian = add_ghost_cells(trian)
+  trian = add_ghost_cells(get_triangulation(f))
   trian_gids = generate_cell_gids(trian)
   cell_to_cellin = _active_aggregates(bgcell_to_bgcellin)
   cell_to_ldofs = cell_ldof_to_mdof(spaces,cell_to_cellin)
@@ -26,7 +24,7 @@ function AgFEMSpace(
 end
 
 function aggregate(strategy,cutgeo::DistributedEmbeddedDiscretization,args...)
-  aggregates,aggregate_owner = distributed_aggregate(strategy,cutgeo,args...)
+  aggregates, aggregate_owner = distributed_aggregate(strategy,cutgeo,args...)
   bgmodel = get_background_model(cutgeo)
   if has_remote_aggregation(bgmodel,aggregates)
     bgmodel = add_remote_aggregates(bgmodel,aggregates,aggregate_owner)
@@ -40,8 +38,8 @@ end
 function distributed_aggregate(
   strategy::AggregateCutCellsByThreshold,
   cut::DistributedEmbeddedDiscretization,
-  in_or_out=IN)
-
+  in_or_out=IN
+)
   geo = get_geometry(cut)
   distributed_aggregate(strategy,cut,geo,in_or_out)
 end
@@ -50,14 +48,13 @@ function distributed_aggregate(
   strategy::AggregateCutCellsByThreshold,
   cut::DistributedEmbeddedDiscretization,
   geo::CSG.Geometry,
-  in_or_out=IN)
-
+  in_or_out=IN
+)
   bgmodel = get_background_model(cut)
   facet_to_inoutcut = compute_bgfacet_to_inoutcut(bgmodel,geo)
   _distributed_aggregate_by_threshold(strategy.threshold,cut,geo,in_or_out,facet_to_inoutcut)
 end
 
-
 function _distributed_aggregate_by_threshold(threshold,cutgeo,geo,loc,facet_to_inoutcut)
   @assert loc in (IN,OUT)
 
@@ -82,15 +79,14 @@ function _distributed_aggregate_by_threshold(threshold,cutgeo,geo,loc,facet_to_i
 
   _distributed_aggregate_by_threshold_barrier(
     threshold,cell_to_unit_cut_meas,facet_to_inoutcut,cell_to_inoutcut,
-    loc,cell_to_coords,cell_to_faces,face_to_cells,gids)
+    loc,cell_to_coords,cell_to_faces,face_to_cells,gids
+  )
 end
 
-
 function _distributed_aggregate_by_threshold_barrier(
   threshold,cell_to_unit_cut_meas,facet_to_inoutcut,cell_to_inoutcut,
-  loc,cell_to_coords,cell_to_faces,face_to_cells,gids)
-
-
+  loc,cell_to_coords,cell_to_faces,face_to_cells,gids
+)
   ocell_to_touched = map(cell_to_unit_cut_meas) do c_to_m
     map(≥,c_to_m,Fill(threshold,length(c_to_m)))
   end
@@ -111,7 +107,6 @@ function _distributed_aggregate_by_threshold_barrier(
   end
 
   cell_to_neig = map(n->zeros(Int32,n),n_cells)
-
   cell_to_root_part = map(collect,local_to_owner(gids))
 
   c1 = map(array_cache,cell_to_faces)
@@ -129,7 +124,8 @@ function _distributed_aggregate_by_threshold_barrier(
       cell_to_faces,
       face_to_cells,
       facet_to_inoutcut,
-      loc)
+      loc
+    )
 
     PVector(cell_to_touched,partition(gids)) |> consistent! |> wait
     PVector(cell_to_neig,partition(gids)) |> consistent! |> wait
@@ -248,8 +244,8 @@ function _find_best_neighbor_from_centroid_distance(
   cell_to_touched,
   cell_to_root_centroid,
   facet_to_inoutcut,
-  loc)
-
+  loc
+)
   faces = getindex!(c1,cell_to_faces,cell)
   dmin = Inf
   T = eltype(eltype(face_to_cells))
@@ -559,6 +555,8 @@ function _local_aggregates(cell_to_gcellin,gcell_to_cell)
   end
 end
 
+# change_bgmodel
+
 function change_bgmodel(cell_to_gcellin,gids::PRange)
   map(change_bgmodel,cell_to_gcellin,local_to_global(gids))
 end
@@ -571,3 +569,87 @@ function change_bgmodel(cell_to_gcellin,ncell_to_gcell)
   end
   ncell_to_gcellin
 end
+
+function change_bgmodel(
+  cutgeo::DistributedEmbeddedDiscretization,
+  model::DistributedDiscreteModel
+)
+  cuts = map(change_bgmodel,local_views(cutgeo),local_views(model))
+  DistributedEmbeddedDiscretization(cuts,model)
+end
+
+function change_bgmodel(
+  cutgeo::DistributedEmbeddedDiscretization,
+  model::DistributedDiscreteModel,
+  cell_to_new_cell
+)
+  cuts = map(change_bgmodel,local_views(cutgeo),local_views(model),cell_to_new_cell)
+  DistributedEmbeddedDiscretization(cuts,model)
+end
+
+function change_bgmodel(
+  cut::EmbeddedDiscretization,
+  newmodel::DiscreteModel,
+  cell_to_newcell=1:num_cells(get_background_model(cut))
+)
+  ls_to_bgc_to_ioc = map(cut.ls_to_bgcell_to_inoutcut) do bgc_to_ioc
+    new_bgc_to_ioc = Vector{Int8}(undef,num_cells(newmodel))
+    new_bgc_to_ioc[cell_to_newcell] = bgc_to_ioc
+    new_bgc_to_ioc
+  end
+  subcells = change_bgmodel(cut.subcells,cell_to_newcell)
+  subfacets = change_bgmodel(cut.subfacets,cell_to_newcell)
+  EmbeddedDiscretization(
+    newmodel,
+    ls_to_bgc_to_ioc,
+    subcells,
+    cut.ls_to_subcell_to_inout,
+    subfacets,
+    cut.ls_to_subfacet_to_inout,
+    cut.oid_to_ls,
+    cut.geo
+  )
+end
+
+function change_bgmodel(
+  cut::EmbeddedFacetDiscretization,
+  newmodel::DiscreteModel,
+  facet_to_newfacet=1:num_facets(get_background_model(cut))
+)
+  nfacets = num_facets(newmodel)
+  ls_to_bgf_to_ioc = map(cut.ls_to_facet_to_inoutcut) do bgf_to_ioc
+    new_bgf_to_ioc = Vector{Int8}(undef,nfacets)
+    new_bgf_to_ioc[facet_to_newfacet] = bgf_to_ioc
+    new_bgf_to_ioc
+  end
+  subfacets = change_bgmodel(cut.subfacets,facet_to_newfacet)
+  EmbeddedFacetDiscretization(
+    newmodel,
+    ls_to_bgf_to_ioc,
+    subfacets,
+    cut.ls_to_subfacet_to_inout,
+    cut.oid_to_ls,
+    cut.geo
+  )
+end
+
+function change_bgmodel(cells::SubCellData,cell_to_newcell)
+  cell_to_bgcell = lazy_map(Reindex(cell_to_newcell),cells.cell_to_bgcell)
+  SubCellData(
+    cells.cell_to_points,
+    collect(Int32,cell_to_bgcell),
+    cells.point_to_coords,
+    cells.point_to_rcoords
+  )
+end
+
+function change_bgmodel(facets::SubFacetData,cell_to_newcell)
+  facet_to_bgcell = lazy_map(Reindex(cell_to_newcell),facets.facet_to_bgcell)
+  SubFacetData(
+    facets.facet_to_points,
+    facets.facet_to_normal,
+    collect(Int32,facet_to_bgcell),
+    facets.point_to_coords,
+    facets.point_to_rcoords
+  )
+end
diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl
index 0e1b602..0f9b55f 100644
--- a/src/Distributed/DistributedDiscreteGeometries.jl
+++ b/src/Distributed/DistributedDiscreteGeometries.jl
@@ -1,44 +1,13 @@
+
 struct DistributedDiscreteGeometry{A} <: GridapType
   geometries::A
 end
 
 local_views(a::DistributedDiscreteGeometry) = a.geometries
 
-# TODO: Is this really necessary? 
-function _get_values_at_owned_coords(φh,model::DistributedDiscreteModel{Dc,Dp}) where {Dc,Dp}
-  @assert DomainStyle(φh) == ReferenceDomain()
-  gids = get_cell_gids(model)
-  values = map(local_views(φh),local_views(model),local_views(gids)) do φh, model, gids
-    own_model = remove_ghost_cells(model,gids)
-    own_cells = get_face_to_parent_face(own_model,Dc)
-
-    trian = get_triangulation(φh)
-    cell_points = get_cell_points(trian)
-    cell_ids = get_cell_node_ids(own_model)
-    cell_values = φh(cell_points)
-    
-    T = eltype(testitem(cell_values))
-    values  = Vector{T}(undef,num_nodes(own_model))
-
-    cell_ids_cache = array_cache(cell_ids)
-    cell_values_cache = array_cache(cell_values)
-    for (ocell,cell) in enumerate(own_cells)
-      ids = getindex!(cell_ids_cache,cell_ids,ocell)
-      vals = getindex!(cell_values_cache,cell_values,cell)
-      values[ids] .= vals
-    end
-    return values
-  end
-  return values
-end
-
 function DiscreteGeometry(φh::CellField,model::DistributedDiscreteModel;name::String="")
-  φ_values = _get_values_at_owned_coords(φh,model)
-  gids = get_cell_gids(model)
-  geometries = map(local_views(model),local_views(gids),local_views(φ_values)) do model,gids,loc_φ
-    ownmodel = remove_ghost_cells(model,gids)
-    point_to_coords = collect1d(get_node_coordinates(ownmodel))
-    DiscreteGeometry(loc_φ,point_to_coords;name)
+  geometries = map(local_views(φh),local_views(model)) do φh, model
+    DiscreteGeometry(φh,model;name)
   end
   DistributedDiscreteGeometry(geometries)
 end
@@ -48,56 +17,45 @@ function distributed_geometry(a::AbstractArray{<:DiscreteGeometry})
 end
 
 function discretize(a::AnalyticalGeometry,model::DistributedDiscreteModel)
-  gids = get_cell_gids(model)
-  geometries = map(local_views(model),local_views(gids)) do model,gids
-    ownmodel = remove_ghost_cells(model,gids)
-    point_to_coords = collect1d(get_node_coordinates(ownmodel))
-    discretize(a,point_to_coords)
+  geometries = map(local_views(model)) do model
+    discretize(a,model)
   end
   DistributedDiscreteGeometry(geometries)
 end
 
-function cut(cutter::Cutter,bgmodel::DistributedDiscreteModel,geom::DistributedDiscreteGeometry)
+function cut(
+  cutter::Cutter,bgmodel::DistributedDiscreteModel,geom::DistributedDiscreteGeometry
+)
   gids = get_cell_gids(bgmodel)
-  cuts = map(local_views(bgmodel),local_views(gids),local_views(geom)) do bgmodel,gids,geom
-    ownmodel = remove_ghost_cells(bgmodel,gids)
-    cutgeo = cut(cutter,ownmodel,geom)
-    change_bgmodel(cutgeo,bgmodel,own_to_local(gids))
+  cuts = map(local_views(bgmodel),local_views(geom)) do bgmodel,geom
+    cut(cutter,bgmodel,geom)
   end
-  consistent_bgcell_to_inoutcut!(cuts,gids)
+  @notimplementedif !isconsistent_bgcell_to_inoutcut(cuts,partition(gids))
   DistributedEmbeddedDiscretization(cuts,bgmodel)
 end
 
-function cut_facets(cutter::Cutter,bgmodel::DistributedDiscreteModel,geom::DistributedDiscreteGeometry)
-  D = map(num_dims,local_views(bgmodel)) |> PartitionedArrays.getany
-  cell_gids = get_cell_gids(bgmodel)
-  facet_gids = get_face_gids(bgmodel,D-1)
-  cuts = map(
-    local_views(bgmodel),
-    local_views(cell_gids),
-    local_views(facet_gids),
-    local_views(geom)) do bgmodel,cell_gids,facet_gids,geom
-    ownmodel = remove_ghost_cells(bgmodel,cell_gids)
-    facet_to_pfacet = get_face_to_parent_face(ownmodel,D-1)
-    cutfacets = cut_facets(cutter,ownmodel,geom)
-    cutfacets = change_bgmodel(cutfacets,bgmodel,facet_to_pfacet)
-    remove_ghost_subfacets(cutfacets,facet_gids)
+function cut_facets(
+  cutter::Cutter,bgmodel::DistributedDiscreteModel{Dc},geom::DistributedDiscreteGeometry
+) where Dc
+  gids = get_face_gids(bgmodel,Dc-1)
+  cuts = map(local_views(bgmodel),local_views(geom)) do bgmodel,geom
+    cut_facets(cutter,bgmodel,geom)
   end
-  consistent_bgfacet_to_inoutcut!(cuts,facet_gids)
+  @notimplementedif !isconsistent_bgcell_to_inoutcut(cuts,partition(gids))
   DistributedEmbeddedDiscretization(cuts,bgmodel)
 end
 
-function distributed_embedded_triangulation(
-  T,
-  cutgeo::DistributedEmbeddedDiscretization,
-  cutinorout,
-  geom::DistributedDiscreteGeometry
-)
-  trians = map(local_views(cutgeo),local_views(geom)) do lcutgeo,lgeom
-    T(lcutgeo,cutinorout,lgeom)
+for TT in (:Triangulation,:SkeletonTriangulation,:BoundaryTriangulation,:EmbeddedBoundary,:GhostSkeleton)
+  @eval begin
+    function $TT(portion,cutgeo::DistributedEmbeddedDiscretization,cutinorout,geom::DistributedDiscreteGeometry)
+      model = get_background_model(cutgeo)
+      gids  = get_cell_gids(model)
+      trians = map(local_views(cutgeo),local_views(geom),partition(gids)) do cutgeo, geom, gids
+        $TT(portion,gids,cutgeo,cutinorout,geom)
+      end
+      DistributedTriangulation(trians,model)
+    end
   end
-  bgmodel = get_background_model(cutgeo)
-  DistributedTriangulation(trians,bgmodel)
 end
 
 function distributed_aggregate(
@@ -111,21 +69,18 @@ function distributed_aggregate(
   _distributed_aggregate_by_threshold(strategy.threshold,cut,geo,in_or_out,facet_to_inoutcut)
 end
 
-function compute_bgcell_to_inoutcut(cutgeo::DistributedEmbeddedDiscretization,geo::DistributedDiscreteGeometry)
-  map(local_views(cutgeo),local_views(geo)) do cutgeo,geo
+function compute_bgcell_to_inoutcut(
+  cutgeo::DistributedEmbeddedDiscretization,geo::DistributedDiscreteGeometry
+)
+  map(local_views(cutgeo),local_views(geo)) do cutgeo, geo
     compute_bgcell_to_inoutcut(cutgeo,geo)
   end
 end
 
 function compute_bgfacet_to_inoutcut(
-  cutter::Cutter,
-  bgmodel::DistributedDiscreteModel,
-  geo::DistributedDiscreteGeometry
+  cutter::Cutter, bgmodel::DistributedDiscreteModel, geo::DistributedDiscreteGeometry
 )
-  gids = get_cell_gids(bgmodel)
-  bgf_to_ioc = map(local_views(bgmodel),local_views(gids),local_views(geo)) do model,gids,geo
-    ownmodel = remove_ghost_cells(model,gids)
-    compute_bgfacet_to_inoutcut(cutter,ownmodel,geo)
+  map(local_views(bgmodel),local_views(geo)) do model, geo
+    compute_bgfacet_to_inoutcut(cutter,model,geo)
   end
-  compute_bgfacet_to_inoutcut(bgmodel,bgf_to_ioc)
 end
\ No newline at end of file
diff --git a/src/Distributed/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl
index e7a8aa1..3fdf7e5 100644
--- a/src/Distributed/DistributedDiscretizations.jl
+++ b/src/Distributed/DistributedDiscretizations.jl
@@ -1,14 +1,14 @@
 
-struct DistributedEmbeddedDiscretization{A,B} <: GridapType
+struct DistributedEmbeddedDiscretization{Dc,Dp,A,B} <: GridapType
   discretizations::A
   model::B
   function DistributedEmbeddedDiscretization(
-    discretizations::AbstractArray{<:AbstractEmbeddedDiscretization},
+    discretizations::AbstractArray{<:AbstractEmbeddedDiscretization{Dc,Dp}},
     model::DistributedDiscreteModel
-  )
+  ) where {Dc,Dp}
     A = typeof(discretizations)
     B = typeof(model)
-    new{A,B}(discretizations,model)
+    new{Dc,Dp,A,B}(discretizations,model)
   end
 end
 
@@ -30,14 +30,12 @@ function cut(bgmodel::DistributedDiscreteModel,args...)
   cut(LevelSetCutter(),bgmodel,args...)
 end
 
-function cut(cutter::Cutter,bgmodel::DistributedDiscreteModel,args...)
-  gids = get_cell_gids(bgmodel)
-  cuts = map(local_views(bgmodel),local_views(gids)) do bgmodel,gids
-    ownmodel = remove_ghost_cells(bgmodel,gids)
-    cutgeo = cut(cutter,ownmodel,args...)
-    change_bgmodel(cutgeo,bgmodel,own_to_local(gids))
+function cut(cutter::Cutter,bgmodel::DistributedDiscreteModel{Dc},args...) where Dc
+  gids = get_face_gids(bgmodel,Dc)
+  cuts = map(local_views(bgmodel)) do bgmodel
+    cut(cutter,bgmodel,args...)
   end
-  consistent_bgcell_to_inoutcut!(cuts,gids)
+  @notimplementedif !isconsistent_bgcell_to_inoutcut(cuts,partition(gids))
   DistributedEmbeddedDiscretization(cuts,bgmodel)
 end
 
@@ -45,306 +43,154 @@ function cut_facets(bgmodel::DistributedDiscreteModel,args...)
   cut_facets(LevelSetCutter(),bgmodel,args...)
 end
 
-function cut_facets(cutter::Cutter,bgmodel::DistributedDiscreteModel,args...)
-  D = map(num_dims,local_views(bgmodel)) |> PartitionedArrays.getany
-  cell_gids = get_cell_gids(bgmodel)
-  facet_gids = get_face_gids(bgmodel,D-1)
-  cuts = map(
-    local_views(bgmodel),
-    local_views(cell_gids),
-    local_views(facet_gids)) do bgmodel,cell_gids,facet_gids
-    ownmodel = remove_ghost_cells(bgmodel,cell_gids)
-    facet_to_pfacet = get_face_to_parent_face(ownmodel,D-1)
-    cutfacets = cut_facets(cutter,ownmodel,args...)
-    cutfacets = change_bgmodel(cutfacets,bgmodel,facet_to_pfacet)
-    remove_ghost_subfacets(cutfacets,facet_gids)
+function cut_facets(cutter::Cutter,bgmodel::DistributedDiscreteModel{Dc},args...) where Dc
+  gids = get_face_gids(bgmodel,Dc-1)
+  cuts = map(local_views(bgmodel)) do bgmodel
+    cut_facets(cutter,bgmodel,args...)
   end
-  consistent_bgfacet_to_inoutcut!(cuts,facet_gids)
+  @notimplementedif !isconsistent_bgcell_to_inoutcut(cuts,partition(gids))
   DistributedEmbeddedDiscretization(cuts,bgmodel)
 end
 
-function Triangulation(
-  cutgeo::DistributedEmbeddedDiscretization,
-  in_or_out::ActiveInOrOut,
-  args...)
-
-  distributed_embedded_triangulation(Triangulation,cutgeo,in_or_out,args...)
-end
-
-function Triangulation(cutgeo::DistributedEmbeddedDiscretization,args...)
-  trian = distributed_embedded_triangulation(Triangulation,cutgeo,args...)
-  remove_ghost_cells(trian)
-end
-
-function EmbeddedBoundary(cutgeo::DistributedEmbeddedDiscretization,args...)
-  trian = distributed_embedded_triangulation(EmbeddedBoundary,cutgeo,args...)
-  remove_ghost_cells(trian)
-end
-
-function SkeletonTriangulation(cutgeo::DistributedEmbeddedDiscretization,args...)
-  trian = distributed_embedded_triangulation(SkeletonTriangulation,cutgeo,args...)
-  remove_ghost_cells(trian)
-end
-
-function BoundaryTriangulation(cutgeo::DistributedEmbeddedDiscretization,args...)
-  trian = distributed_embedded_triangulation(BoundaryTriangulation,cutgeo,args...)
-  remove_ghost_cells(trian)
-end
-
-function distributed_embedded_triangulation(
-  T,
-  cutgeo::DistributedEmbeddedDiscretization,
-  args...)
+# Note on distributed triangulations:
+#
+# - We allow for more one argument, `portion`, which allows the user to filter
+# some of the cells/faces. In particular, this is used to remove ghosts from the
+# local triangulations.
+# - The default for `portion` is `NoGhost()`, wich filters out all ghost cells, except
+# when we have the argument `in_or_out`.
 
-  trians = map(local_views(cutgeo)) do lcutgeo
-    T(lcutgeo,args...)
-  end
-  bgmodel = get_background_model(cutgeo)
-  DistributedTriangulation(trians,bgmodel)
-end
-
-function compute_bgfacet_to_inoutcut(
-  bgmodel::DistributedDiscreteModel,
-  bgf_to_ioc::AbstractArray{<:AbstractVector})
-
-  D = num_dims(eltype(local_views(bgmodel)))
-  gids = get_cell_gids(bgmodel)
-  bgf_to_ioc = map(
-    local_views(bgmodel),
-    local_views(gids),
-    bgf_to_ioc) do bgmodel,gids,bgf_to_ioc
-
-    ownmodel = remove_ghost_cells(bgmodel,gids)
-    f_to_pf = Gridap.Geometry.get_face_to_parent_face(ownmodel,D-1)
-    _bgf_to_ioc = Vector{eltype(bgf_to_ioc)}(undef,num_faces(bgmodel,D-1))
-    _bgf_to_ioc[f_to_pf] .= bgf_to_ioc
-    _bgf_to_ioc
-  end
-  facet_gids = get_face_gids(bgmodel,D-1)
-  pbgf_to_ioc = PVector(bgf_to_ioc,partition(facet_gids))
-  consistent!(pbgf_to_ioc) |> wait
-  local_values(pbgf_to_ioc)
-end
-
-function compute_bgfacet_to_inoutcut(
-  cutter::Cutter,
-  bgmodel::DistributedDiscreteModel,
-  geo
+function Triangulation(
+  cutgeo::DistributedEmbeddedDiscretization,in_or_out::ActiveInOrOut,args...
 )
-  gids = get_cell_gids(bgmodel)
-  bgf_to_ioc = map(local_views(bgmodel),local_views(gids)) do model,gids
-    ownmodel = remove_ghost_cells(model,gids)
-    compute_bgfacet_to_inoutcut(cutter,ownmodel,geo)
-  end
-  compute_bgfacet_to_inoutcut(bgmodel,bgf_to_ioc)
+  Triangulation(WithGhost(),cutgeo,in_or_out,args...)
 end
 
-function compute_bgfacet_to_inoutcut(bgmodel::DistributedDiscreteModel,args...)
-  cutter = LevelSetCutter()
-  compute_bgfacet_to_inoutcut(cutter,bgmodel,args...)
-end
-
-function compute_bgcell_to_inoutcut(cutgeo::DistributedEmbeddedDiscretization,args...)
-  map(local_views(cutgeo)) do cutgeo
-    compute_bgcell_to_inoutcut(cutgeo,args...)
-  end
-end
+for TT in (:Triangulation,:SkeletonTriangulation,:BoundaryTriangulation,:EmbeddedBoundary,:GhostSkeleton)
+  @eval begin
+    function $TT(cutgeo::DistributedEmbeddedDiscretization,args...)
+      $TT(NoGhost(),cutgeo,args...)
+    end
 
-function compute_bgfacet_to_inoutcut(cutgeo::DistributedEmbeddedDiscretization,args...)
-  map(local_views(cutgeo)) do cutgeo
-    compute_bgfacet_to_inoutcut(cutgeo,args...)
-  end
-end
+    function $TT(portion,cutgeo::DistributedEmbeddedDiscretization,args...)
+      model = get_background_model(cutgeo)
+      gids  = get_cell_gids(model)
+      trians = map(local_views(cutgeo),partition(gids)) do cutgeo, gids
+        $TT(portion,gids,cutgeo,args...)
+      end
+      DistributedTriangulation(trians,model)
+    end
 
-function remove_ghost_cells(trian::DistributedTriangulation)
-  model = get_background_model(trian)
-  gids = get_cell_gids(model)
-  trians = map(local_views(trian),local_views(gids)) do trian,gids
-    remove_ghost_cells(trian,gids)
+    function $TT(portion,gids::AbstractLocalIndices,cutgeo::AbstractEmbeddedDiscretization,args...)
+      trian = $TT(cutgeo,args...)
+      filter_cells_when_needed(portion,gids,trian)
+    end
   end
-  DistributedTriangulation(trians,model)
 end
 
+# TODO: This should go to GridapDistributed
 function remove_ghost_cells(trian::AppendedTriangulation,gids)
   a = remove_ghost_cells(trian.a,gids)
   b = remove_ghost_cells(trian.b,gids)
   lazy_append(a,b)
 end
 
-function remove_ghost_cells(trian::SubFacetTriangulation,gids)
-  model = get_background_model(trian)
-  D     = num_cell_dims(model)
-  glue  = get_glue(trian,Val{D}())
+function remove_ghost_cells(trian::SubFacetTriangulation{Df,Dc},gids) where {Df,Dc}
+  glue  = get_glue(trian,Val{Dc}())
   remove_ghost_cells(glue,trian,gids)
 end
 
-function remove_ghost_cells(model::DiscreteModel,gids::AbstractLocalIndices)
-  DiscreteModelPortion(model,own_to_local(gids))
-end
-
-function consistent_bgcell_to_inoutcut!(
-  cuts::AbstractArray{<:AbstractEmbeddedDiscretization},
-  gids::PRange)
-
-  ls_to_bgcell_to_inoutcut = map(get_ls_to_bgcell_to_inoutcut,cuts)
-  _consistent!(ls_to_bgcell_to_inoutcut,gids)
-end
-
-function get_ls_to_bgcell_to_inoutcut(cut::EmbeddedDiscretization)
-  cut.ls_to_bgcell_to_inoutcut
-end
-
-function consistent_bgfacet_to_inoutcut!(
-  cuts::AbstractArray{<:AbstractEmbeddedDiscretization},
-  gids::PRange)
-
-  ls_to_bgfacet_to_inoutcut = map(get_ls_to_bgfacet_to_inoutcut,cuts)
-  _consistent!(ls_to_bgfacet_to_inoutcut,gids)
-end
-
-function get_ls_to_bgfacet_to_inoutcut(cut::EmbeddedFacetDiscretization)
-  cut.ls_to_facet_to_inoutcut
-end
-
-function _consistent!(
-  p_to_i_to_a::AbstractArray{<:Vector{<:Vector}},
-  prange::PRange)
-
-  n = map(length,p_to_i_to_a) |> PartitionedArrays.getany
-  for i in 1:n
-    p_to_a = map(i_to_a->i_to_a[i],p_to_i_to_a)
-    PVector(p_to_a,partition(prange)) |> consistent! |> wait
-    map(p_to_a,p_to_i_to_a) do p_to_a,p_to_ia
-      copyto!(p_to_ia[i],p_to_a)
-    end
+function remove_ghost_subfacets(cut::EmbeddedFacetDiscretization,facet_gids)
+  bgfacet_mask = map(!iszero,local_to_owner(facet_gids))
+  subfacet_mask = map(Reindex(bgfacet_mask),cut.subfacets.cell_to_bgcell)
+  new_subfacets = findall(subfacet_mask)
+  subfacets = SubCellData(cut.subfacets,new_subfacets)
+  ls_to_subfacet_to_inout = map(cut.ls_to_subfacet_to_inout) do sf_to_io
+    map(Reindex(sf_to_io),new_subfacets)
   end
+  EmbeddedFacetDiscretization(
+    cut.bgmodel,
+    cut.ls_to_facet_to_inoutcut,
+    subfacets,
+    ls_to_subfacet_to_inout,
+    cut.oid_to_ls,
+    cut.geo
+  )
 end
 
-function change_bgmodel(
-  cutgeo::DistributedEmbeddedDiscretization,
-  model::DistributedDiscreteModel,
-  args...)
-
-  cuts = _change_bgmodels(cutgeo,model,args...)
-  gids = get_cell_gids(model)
-  ls_to_bgcell_to_inoutcut = map(c->c.ls_to_bgcell_to_inoutcut,cuts)
-  _consistent!(ls_to_bgcell_to_inoutcut,gids)
-  DistributedEmbeddedDiscretization(cuts,model)
-end
-
-function change_bgmodel(
-  cutgeo::DistributedEmbeddedDiscretization{<:AbstractArray{<:EmbeddedFacetDiscretization}},
-  model::DistributedDiscreteModel,
-  args...)
+# Distributed InOutCut flag methods
 
-  D = map(num_dims,local_views(model)) |> PartitionedArrays.getany
-  cuts = _change_bgmodels(cutgeo,model,args...)
-  gids = get_face_gids(model,D-1)
-  ls_to_facet_to_inoutcut = map(c->c.ls_to_facet_to_inoutcut,cuts)
-  _consistent!(ls_to_facet_to_inoutcut,gids)
-  DistributedEmbeddedDiscretization(cuts,model)
+#     isconsistent_bgcell_to_inoutcut(cut::DistributedEmbeddedDiscretization)
+#     isconsistent_bgcell_to_inoutcut(cuts::AbstractArray{<:AbstractEmbeddedDiscretization},indices)
+#
+# Returns true if the local `ls_to_bgcell_to_inoutcut` arrays are consistent
+# accross processors.
+function isconsistent_bgcell_to_inoutcut(
+  cut::DistributedEmbeddedDiscretization{Dc}
+) where Dc
+  model = get_background_model(cut)
+  gids = get_face_gids(model,Dc)
+  isconsistent_bgcell_to_inoutcut(local_views(cut),partition(gids))
 end
 
-function _change_bgmodels(
-  cutgeo::DistributedEmbeddedDiscretization,
-  model::DistributedDiscreteModel,
-  cell_to_newcell
+function isconsistent_bgcell_to_inoutcut(
+  cuts::AbstractArray{<:AbstractEmbeddedDiscretization},indices::AbstractArray
 )
-  map(local_views(cutgeo),local_views(model),cell_to_newcell) do c,m,c_to_nc
-    change_bgmodel(c,m,c_to_nc)
+  get_inoutcut(cut::EmbeddedDiscretization) = Tuple(cut.ls_to_bgcell_to_inoutcut)
+  get_inoutcut(cut::EmbeddedFacetDiscretization) = Tuple(cut.ls_to_facet_to_inoutcut)
+  ls_to_bgcell_to_inoutcut = tuple_of_arrays(map(get_inoutcut,cuts))
+  return isconsistent_bgcell_to_inoutcut(ls_to_bgcell_to_inoutcut,indices)
+end
+
+function isconsistent_bgcell_to_inoutcut(
+  ls_to_bgcell_to_inoutcut::NTuple{N,<:AbstractArray{<:Vector}},indices::AbstractArray
+) where N
+  for bgcell_to_inoutcut in ls_to_bgcell_to_inoutcut
+    if !isconsistent_bgcell_to_inoutcut(bgcell_to_inoutcut,indices)
+      return false
+    end
   end
+  return true
 end
 
-function _change_bgmodels(
-  cutgeo::DistributedEmbeddedDiscretization,
-  model::DistributedDiscreteModel
+function isconsistent_bgcell_to_inoutcut(
+  bgcell_to_inoutcut::AbstractArray{<:Vector},indices::AbstractArray
 )
-  map(local_views(cutgeo),local_views(model)) do c,m
-    change_bgmodel(c,m)
+  # TODO: Some allocations can be avoided by going to the low-level communication API
+  ref = map(copy,bgcell_to_inoutcut)
+  wait(consistent!(PVector(ref,indices)))
+  is_consistent = map(bgcell_to_inoutcut,ref) do bgcell_to_inoutcut,ref
+    bgcell_to_inoutcut == ref
   end
+  return reduce(&,is_consistent,init=true)
 end
 
-function change_bgmodel(
-  cut::EmbeddedDiscretization,
-  newmodel::DiscreteModel,
-  cell_to_newcell=1:num_cells(get_background_model(cut))
-)
-  ls_to_bgc_to_ioc = map(cut.ls_to_bgcell_to_inoutcut) do bgc_to_ioc
-    new_bgc_to_ioc = Vector{Int8}(undef,num_cells(newmodel))
-    new_bgc_to_ioc[cell_to_newcell] = bgc_to_ioc
-    new_bgc_to_ioc
-  end
-  subcells = change_bgmodel(cut.subcells,cell_to_newcell)
-  subfacets = change_bgmodel(cut.subfacets,cell_to_newcell)
-  EmbeddedDiscretization(
-    newmodel,
-    ls_to_bgc_to_ioc,
-    subcells,
-    cut.ls_to_subcell_to_inout,
-    subfacets,
-    cut.ls_to_subfacet_to_inout,
-    cut.oid_to_ls,
-    cut.geo)
+# TODO: Should we check for consistency here?
+function compute_bgfacet_to_inoutcut(bgmodel::DistributedDiscreteModel,args...)
+  cutter = LevelSetCutter()
+  compute_bgfacet_to_inoutcut(cutter,bgmodel,args...)
 end
 
-function change_bgmodel(
-  cut::EmbeddedFacetDiscretization,
-  newmodel::DiscreteModel,
-  facet_to_newfacet=1:num_facets(get_background_model(cut))
-)
-  nfacets = num_facets(newmodel)
-  ls_to_bgf_to_ioc = map(cut.ls_to_facet_to_inoutcut) do bgf_to_ioc
-    new_bgf_to_ioc = Vector{Int8}(undef,nfacets)
-    new_bgf_to_ioc[facet_to_newfacet] = bgf_to_ioc
-    new_bgf_to_ioc
+function compute_bgfacet_to_inoutcut(cutter::Cutter,bgmodel::DistributedDiscreteModel,args...)
+  map(local_views(bgmodel)) do bgmodel
+    compute_bgfacet_to_inoutcut(cutter,bgmodel,args...)
   end
-  subfacets = change_bgmodel(cut.subfacets,facet_to_newfacet)
-  EmbeddedFacetDiscretization(
-    newmodel,
-    ls_to_bgf_to_ioc,
-    subfacets,
-    cut.ls_to_subfacet_to_inout,
-    cut.oid_to_ls,
-    cut.geo
-  )
-end
-
-function change_bgmodel(cells::SubCellData,cell_to_newcell)
-  cell_to_bgcell = lazy_map(Reindex(cell_to_newcell),cells.cell_to_bgcell)
-  SubCellData(
-    cells.cell_to_points,
-    collect(Int32,cell_to_bgcell),
-    cells.point_to_coords,
-    cells.point_to_rcoords)
 end
 
-function change_bgmodel(facets::SubFacetData,cell_to_newcell)
-  facet_to_bgcell = lazy_map(Reindex(cell_to_newcell),facets.facet_to_bgcell)
-  SubFacetData(
-    facets.facet_to_points,
-    facets.facet_to_normal,
-    collect(Int32,facet_to_bgcell),
-    facets.point_to_coords,
-    facets.point_to_rcoords)
+function compute_bgcell_to_inoutcut(cutgeo::DistributedEmbeddedDiscretization,args...)
+  map(local_views(cutgeo)) do cutgeo
+    compute_bgcell_to_inoutcut(cutgeo,args...)
+  end
 end
 
-function remove_ghost_subfacets(cut::EmbeddedFacetDiscretization,facet_gids)
-  bgfacet_mask = map(!iszero,local_to_owner(facet_gids))
-  subfacet_mask = map(Reindex(bgfacet_mask),cut.subfacets.cell_to_bgcell)
-  new_subfacets = findall(subfacet_mask)
-  subfacets = SubCellData(cut.subfacets,new_subfacets)
-  ls_to_subfacet_to_inout = map(cut.ls_to_subfacet_to_inout) do sf_to_io
-    map(Reindex(sf_to_io),new_subfacets)
+function compute_bgfacet_to_inoutcut(cutgeo::DistributedEmbeddedDiscretization,args...)
+  map(local_views(cutgeo)) do cutgeo
+    compute_bgfacet_to_inoutcut(cutgeo,args...)
   end
-  EmbeddedFacetDiscretization(
-    cut.bgmodel,
-    cut.ls_to_facet_to_inoutcut,
-    subfacets,
-    ls_to_subfacet_to_inout,
-    cut.oid_to_ls,
-    cut.geo)
 end
 
+# AMR
+
 function compute_redistribute_wights(
   cut::DistributedEmbeddedDiscretization,
   args...)
diff --git a/src/Distributed/DistributedQuadratures.jl b/src/Distributed/DistributedQuadratures.jl
index 06a40a4..634b778 100644
--- a/src/Distributed/DistributedQuadratures.jl
+++ b/src/Distributed/DistributedQuadratures.jl
@@ -1,18 +1,18 @@
 
 function CellData.Measure(
-  t::DistributedTriangulation,
+  trian::DistributedTriangulation,
   quad::Tuple{MomentFitted,Vararg};
   kwargs...)
   @notimplemented
   name, _args, _kwargs = quad
   cut,cutfacets,_args... = _args
-  t = remove_ghost_cells(t)
-  measures = map(
-    local_views(t),
-    local_views(cut),
-    local_views(cutfacets)) do trian,cut,cutfacets
-      quad = name, (cut,cutfacets,_args...), _kwargs
-      Measure(trian,quad;kwargs...)
+
+  model = get_background_model(trian)
+  gids  = get_cell_gids(model)
+  trian = remove_ghost_cells(trian,gids)
+  measures = map(local_views(trian),local_views(cut),local_views(cutfacets)) do trian,cut,cutfacets
+    quad = name, (cut,cutfacets,_args...), _kwargs
+    Measure(trian,quad;kwargs...)
   end
   DistributedMeasure(measures)
 end
diff --git a/src/Distributed/DistributedSubFacetTriangulations.jl b/src/Distributed/DistributedSubFacetTriangulations.jl
new file mode 100644
index 0000000..d47c350
--- /dev/null
+++ b/src/Distributed/DistributedSubFacetTriangulations.jl
@@ -0,0 +1,63 @@
+
+const DistributedSubFacetTriangulation{Df,Dc} = DistributedTriangulation{Df,Dc,<:AbstractArray{<:Union{SubFacetTriangulation{Df,Dc},TriangulationView{Df,Dc,<:SubFacetTriangulation{Df,Dc}}}}}
+
+# Each cut facet belongs to the background cell containing it. So we can generate 
+# ownership information for the cut facets from the background cell gids.
+function GridapDistributed.generate_cell_gids(
+  trian::DistributedSubFacetTriangulation{Df,Dc},
+) where {Df,Dc}
+  model = get_background_model(trian)
+  cgids = get_cell_gids(model)
+  
+  n_lfacets = map(num_cells,local_views(trian))
+  first_gid = scan(+,n_lfacets,type=:exclusive,init=one(eltype(n_lfacets)))
+  n_facets = reduce(+,n_lfacets,init=zero(eltype(n_lfacets)))
+
+  fgids = map(local_views(trian),partition(cgids),first_gid,n_lfacets) do trian, cgids, first_gid, n_lfacets
+    glue = get_glue(trian,Val(Dc)) # Glue from cut facets to background cells 
+    facet_to_bgcell = glue.tface_to_mface
+
+    facet_to_gid = collect(first_gid:(first_gid+n_lfacets-1))
+    facet_to_owner = local_to_owner(cgids)[facet_to_bgcell]
+    LocalIndices(n_facets,part_id(cgids),facet_to_gid,facet_to_owner)
+  end
+  return PRange(fgids)
+end
+
+function GridapDistributed.add_ghost_cells(
+  trian::DistributedSubFacetTriangulation{Df,Dc},
+) where {Df,Dc}
+
+  # In this case, we already have all ghost facets
+  if eltype(local_views(trian)) <: SubFacetTriangulation
+    return trian
+  end
+
+  # First, we create a new Triangulation containing all the cut facets
+  model = get_background_model(trian)
+  bgtrians, facet_to_bgfacet = map(local_views(trian)) do trian
+    @assert isa(trian,TriangulationView)
+    trian.parent, trian.cell_to_parent_cell
+  end |> tuple_of_arrays
+  bgtrian = DistributedTriangulation(bgtrians,model)
+  fgids = partition(generate_cell_gids(bgtrian))
+
+  # Exchange info about cut facets
+  inside_facets = map(fgids,facet_to_bgfacet) do fgids, facet_to_bgfacet
+    inside_facets = falses(local_length(fgids))
+    inside_facets[facet_to_bgfacet] .= true
+    return inside_facets
+  end
+  wait(consistent!(PVector(inside_facets,fgids))) # Exchange information
+
+  # Return ghosted Triangulation
+  covers_all = reduce(&,map(all,inside_facets),init=true)
+  if covers_all
+    ghosted_trian = bgtrian
+  else
+    ghosted_trian = DistributedTriangulation(
+      map(TriangulationView,bgtrians,inside_facets), model
+    )
+  end
+  return ghosted_trian
+end
diff --git a/src/Interfaces/EmbeddedDiscretizations.jl b/src/Interfaces/EmbeddedDiscretizations.jl
index 9dcf623..61d2e75 100644
--- a/src/Interfaces/EmbeddedDiscretizations.jl
+++ b/src/Interfaces/EmbeddedDiscretizations.jl
@@ -1,12 +1,12 @@
 
-abstract type AbstractEmbeddedDiscretization <: GridapType end
+abstract type AbstractEmbeddedDiscretization{Dc,Dp} <: GridapType end
 
-struct EmbeddedDiscretization{Dp,T} <: AbstractEmbeddedDiscretization
+struct EmbeddedDiscretization{Dc,T} <: AbstractEmbeddedDiscretization{Dc,Dc}
   bgmodel::DiscreteModel
   ls_to_bgcell_to_inoutcut::Vector{Vector{Int8}}
-  subcells::SubCellData{Dp,Dp,T}
+  subcells::SubCellData{Dc,Dc,T}
   ls_to_subcell_to_inout::Vector{Vector{Int8}}
-  subfacets::SubFacetData{Dp,T}
+  subfacets::SubFacetData{Dc,T}
   ls_to_subfacet_to_inout::Vector{Vector{Int8}}
   oid_to_ls::Dict{UInt,Int}
   geo::CSG.Geometry
diff --git a/src/Interfaces/EmbeddedFacetDiscretizations.jl b/src/Interfaces/EmbeddedFacetDiscretizations.jl
index 4086fb6..b521906 100644
--- a/src/Interfaces/EmbeddedFacetDiscretizations.jl
+++ b/src/Interfaces/EmbeddedFacetDiscretizations.jl
@@ -1,5 +1,5 @@
 
-struct EmbeddedFacetDiscretization{Dc,Dp,T} <: AbstractEmbeddedDiscretization
+struct EmbeddedFacetDiscretization{Dc,Dp,T} <: AbstractEmbeddedDiscretization{Dc,Dp}
   bgmodel::DiscreteModel{Dp,Dp}
   ls_to_facet_to_inoutcut::Vector{Vector{Int8}}
   subfacets::SubCellData{Dc,Dp,T}