From 3185beb8c5e2af57af4bc9009d9be659fdd455b3 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 14 Nov 2024 14:50:38 +1100
Subject: [PATCH 01/11] Started refactoring the distributed code

---
 src/Distributed/Distributed.jl                |  13 +-
 .../DistributedDiscreteGeometries.jl          | 114 ++----
 src/Distributed/DistributedDiscretizations.jl | 375 ++++++------------
 src/Interfaces/EmbeddedDiscretizations.jl     |   8 +-
 .../EmbeddedFacetDiscretizations.jl           |   2 +-
 5 files changed, 156 insertions(+), 356 deletions(-)

diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl
index ee2ea62d..e204227b 100644
--- a/src/Distributed/Distributed.jl
+++ b/src/Distributed/Distributed.jl
@@ -31,14 +31,11 @@ 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
diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl
index 0e1b6021..83f96ad9 100644
--- a/src/Distributed/DistributedDiscreteGeometries.jl
+++ b/src/Distributed/DistributedDiscreteGeometries.jl
@@ -4,41 +4,9 @@ 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 +16,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)
+  @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 +68,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 e7a8aa1b..3c84be45 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,124 +43,51 @@ 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: 
 
-  trians = map(local_views(cutgeo)) do lcutgeo
-    T(lcutgeo,args...)
-  end
-  bgmodel = get_background_model(cutgeo)
-  DistributedTriangulation(trians,bgmodel)
-end
+- 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`. 
+"""
 
-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)
+  @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
 
 function remove_ghost_cells(trian::AppendedTriangulation,gids)
@@ -171,180 +96,104 @@ function remove_ghost_cells(trian::AppendedTriangulation,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
+# Distributed InOutCut flag methods
 
-function change_bgmodel(
-  cutgeo::DistributedEmbeddedDiscretization{<:AbstractArray{<:EmbeddedFacetDiscretization}},
-  model::DistributedDiscreteModel,
-  args...)
+"""
+    isconsistent_bgcell_to_inoutcut(cut::DistributedEmbeddedDiscretization)
+    isconsistent_bgcell_to_inoutcut(cuts::AbstractArray{<:AbstractEmbeddedDiscretization},indices)
 
-  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)
+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/Interfaces/EmbeddedDiscretizations.jl b/src/Interfaces/EmbeddedDiscretizations.jl
index 9dcf6237..61d2e759 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 4086fb69..b521906c 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}

From 1d82808cd65bcff9557a1d3d57ef071db992f3dc Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 14 Nov 2024 15:38:58 +1100
Subject: [PATCH 02/11] Added add_ghost_cells and generate_cell_gids for
 distributed SubFacetTriangulations

---
 src/Distributed/Distributed.jl                |  2 +
 .../DistributedDiscreteGeometries.jl          |  1 +
 src/Distributed/DistributedDiscretizations.jl | 17 +++--
 .../DistributedSubFacetTriangulations.jl      | 63 +++++++++++++++++++
 4 files changed, 74 insertions(+), 9 deletions(-)
 create mode 100644 src/Distributed/DistributedSubFacetTriangulations.jl

diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl
index e204227b..6d881bba 100644
--- a/src/Distributed/Distributed.jl
+++ b/src/Distributed/Distributed.jl
@@ -58,6 +58,8 @@ include("DistributedDiscretizations.jl")
 
 include("DistributedDiscreteGeometries.jl")
 
+include("DistributedSubFacetTriangulations.jl")
+
 include("DistributedAgFEM.jl")
 
 include("DistributedQuadratures.jl")
diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl
index 83f96ad9..f196c6e9 100644
--- a/src/Distributed/DistributedDiscreteGeometries.jl
+++ b/src/Distributed/DistributedDiscreteGeometries.jl
@@ -1,3 +1,4 @@
+
 struct DistributedDiscreteGeometry{A} <: GridapType
   geometries::A
 end
diff --git a/src/Distributed/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl
index 3c84be45..a35de042 100644
--- a/src/Distributed/DistributedDiscretizations.jl
+++ b/src/Distributed/DistributedDiscretizations.jl
@@ -52,15 +52,13 @@ function cut_facets(cutter::Cutter,bgmodel::DistributedDiscreteModel{Dc},args...
   DistributedEmbeddedDiscretization(cuts,bgmodel)
 end
 
-"""
-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`. 
-"""
+# 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`.
 
 function Triangulation(
   cutgeo::DistributedEmbeddedDiscretization,in_or_out::ActiveInOrOut,args...
@@ -90,6 +88,7 @@ for TT in (:Triangulation,:SkeletonTriangulation,:BoundaryTriangulation,:Embedde
   end
 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)
diff --git a/src/Distributed/DistributedSubFacetTriangulations.jl b/src/Distributed/DistributedSubFacetTriangulations.jl
new file mode 100644
index 00000000..d47c3501
--- /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

From 946cdc120bf72ce251e1d608cd09abdbcf3a9d66 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 14 Nov 2024 15:45:49 +1100
Subject: [PATCH 03/11] Fixed distributed moment-fitted quads (maybe?)

---
 src/Distributed/DistributedQuadratures.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/Distributed/DistributedQuadratures.jl b/src/Distributed/DistributedQuadratures.jl
index 06a40a4c..634b7786 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

From 7eac5af045a120f5b10a6dde3132b2f00d282d62 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 14 Nov 2024 16:38:16 +1100
Subject: [PATCH 04/11] Fixes

---
 src/Distributed/DistributedAgFEM.jl           | 122 +++++++++++++++---
 src/Distributed/DistributedDiscretizations.jl |  12 +-
 2 files changed, 107 insertions(+), 27 deletions(-)

diff --git a/src/Distributed/DistributedAgFEM.jl b/src/Distributed/DistributedAgFEM.jl
index 4501ba76..cdb65ab1 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/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl
index a35de042..5d6a22ba 100644
--- a/src/Distributed/DistributedDiscretizations.jl
+++ b/src/Distributed/DistributedDiscretizations.jl
@@ -120,13 +120,11 @@ end
 
 # Distributed InOutCut flag methods
 
-"""
-    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. 
-"""
+#     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

From af3720f38dbee38a3d20939648f1660dca431d43 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 14 Nov 2024 16:44:42 +1100
Subject: [PATCH 05/11] Minor

---
 src/Distributed/Distributed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl
index 6d881bba..47d92880 100644
--- a/src/Distributed/Distributed.jl
+++ b/src/Distributed/Distributed.jl
@@ -26,7 +26,7 @@ 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

From 392b451f8f72d551bd4f665ea40941eacf966f82 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 14 Nov 2024 21:55:26 +1100
Subject: [PATCH 06/11] Updated NEWS

---
 NEWS.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/NEWS.md b/NEWS.md
index 7d7bb3a4..325955f8 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
 

From 6850b686cc41cf5a018436377fc35efa63e60ed2 Mon Sep 17 00:00:00 2001
From: zjwegert <60646897+zjwegert@users.noreply.github.com>
Date: Mon, 25 Nov 2024 13:19:20 +1000
Subject: [PATCH 07/11] Distributed ghost skeleton

---
 src/Distributed/Distributed.jl                |  1 +
 .../DistributedDiscreteGeometries.jl          |  2 +-
 src/Distributed/DistributedDiscretizations.jl | 22 +++++++++----------
 3 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl
index 47d92880..b4cf0b8e 100644
--- a/src/Distributed/Distributed.jl
+++ b/src/Distributed/Distributed.jl
@@ -45,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
diff --git a/src/Distributed/DistributedDiscreteGeometries.jl b/src/Distributed/DistributedDiscreteGeometries.jl
index f196c6e9..0f9b55f2 100644
--- a/src/Distributed/DistributedDiscreteGeometries.jl
+++ b/src/Distributed/DistributedDiscreteGeometries.jl
@@ -45,7 +45,7 @@ function cut_facets(
   DistributedEmbeddedDiscretization(cuts,bgmodel)
 end
 
-for TT in (:Triangulation,:SkeletonTriangulation,:BoundaryTriangulation,:EmbeddedBoundary)
+for TT in (:Triangulation,:SkeletonTriangulation,:BoundaryTriangulation,:EmbeddedBoundary,:GhostSkeleton)
   @eval begin
     function $TT(portion,cutgeo::DistributedEmbeddedDiscretization,cutinorout,geom::DistributedDiscreteGeometry)
       model = get_background_model(cutgeo)
diff --git a/src/Distributed/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl
index 5d6a22ba..3fdf7e5d 100644
--- a/src/Distributed/DistributedDiscretizations.jl
+++ b/src/Distributed/DistributedDiscretizations.jl
@@ -52,12 +52,12 @@ function cut_facets(cutter::Cutter,bgmodel::DistributedDiscreteModel{Dc},args...
   DistributedEmbeddedDiscretization(cuts,bgmodel)
 end
 
-# 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 
+# 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`.
 
 function Triangulation(
@@ -66,7 +66,7 @@ function Triangulation(
   Triangulation(WithGhost(),cutgeo,in_or_out,args...)
 end
 
-for TT in (:Triangulation,:SkeletonTriangulation,:BoundaryTriangulation,:EmbeddedBoundary)
+for TT in (:Triangulation,:SkeletonTriangulation,:BoundaryTriangulation,:EmbeddedBoundary,:GhostSkeleton)
   @eval begin
     function $TT(cutgeo::DistributedEmbeddedDiscretization,args...)
       $TT(NoGhost(),cutgeo,args...)
@@ -122,9 +122,9 @@ end
 
 #     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. 
+#
+# 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
@@ -189,7 +189,7 @@ function compute_bgfacet_to_inoutcut(cutgeo::DistributedEmbeddedDiscretization,a
   end
 end
 
-# AMR 
+# AMR
 
 function compute_redistribute_wights(
   cut::DistributedEmbeddedDiscretization,

From a8bb48d051be7f0b49f1e4da09cc6d1b817a3435 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Thu, 27 Feb 2025 17:09:49 +1100
Subject: [PATCH 08/11] Minor

---
 src/Distributed/DistributedDiscretizations.jl |  8 ++++----
 src/Interfaces/EmbeddedDiscretizations.jl     |  4 ++--
 .../EmbeddedFacetDiscretizations.jl           |  2 +-
 src/LevelSetCutters/CutTriangulations.jl      | 20 +++++++++----------
 4 files changed, 16 insertions(+), 18 deletions(-)

diff --git a/src/Distributed/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl
index 5d6a22ba..0a43c2f6 100644
--- a/src/Distributed/DistributedDiscretizations.jl
+++ b/src/Distributed/DistributedDiscretizations.jl
@@ -1,14 +1,14 @@
 
-struct DistributedEmbeddedDiscretization{Dc,Dp,A,B} <: GridapType
+struct DistributedEmbeddedDiscretization{A,B} <: GridapType
   discretizations::A
   model::B
   function DistributedEmbeddedDiscretization(
-    discretizations::AbstractArray{<:AbstractEmbeddedDiscretization{Dc,Dp}},
+    discretizations::AbstractArray{<:AbstractEmbeddedDiscretization},
     model::DistributedDiscreteModel
-  ) where {Dc,Dp}
+  )
     A = typeof(discretizations)
     B = typeof(model)
-    new{Dc,Dp,A,B}(discretizations,model)
+    new{A,B}(discretizations,model)
   end
 end
 
diff --git a/src/Interfaces/EmbeddedDiscretizations.jl b/src/Interfaces/EmbeddedDiscretizations.jl
index 61d2e759..5a6a14b3 100644
--- a/src/Interfaces/EmbeddedDiscretizations.jl
+++ b/src/Interfaces/EmbeddedDiscretizations.jl
@@ -1,7 +1,7 @@
 
-abstract type AbstractEmbeddedDiscretization{Dc,Dp} <: GridapType end
+abstract type AbstractEmbeddedDiscretization <: GridapType end
 
-struct EmbeddedDiscretization{Dc,T} <: AbstractEmbeddedDiscretization{Dc,Dc}
+struct EmbeddedDiscretization{Dc,T} <: AbstractEmbeddedDiscretization
   bgmodel::DiscreteModel
   ls_to_bgcell_to_inoutcut::Vector{Vector{Int8}}
   subcells::SubCellData{Dc,Dc,T}
diff --git a/src/Interfaces/EmbeddedFacetDiscretizations.jl b/src/Interfaces/EmbeddedFacetDiscretizations.jl
index b521906c..4086fb69 100644
--- a/src/Interfaces/EmbeddedFacetDiscretizations.jl
+++ b/src/Interfaces/EmbeddedFacetDiscretizations.jl
@@ -1,5 +1,5 @@
 
-struct EmbeddedFacetDiscretization{Dc,Dp,T} <: AbstractEmbeddedDiscretization{Dc,Dp}
+struct EmbeddedFacetDiscretization{Dc,Dp,T} <: AbstractEmbeddedDiscretization
   bgmodel::DiscreteModel{Dp,Dp}
   ls_to_facet_to_inoutcut::Vector{Vector{Int8}}
   subfacets::SubCellData{Dc,Dp,T}
diff --git a/src/LevelSetCutters/CutTriangulations.jl b/src/LevelSetCutters/CutTriangulations.jl
index f1bd1f69..17068706 100644
--- a/src/LevelSetCutters/CutTriangulations.jl
+++ b/src/LevelSetCutters/CutTriangulations.jl
@@ -149,10 +149,14 @@ function cut_sub_triangulation(m::CutTriangulation, mpoint_to_value)
 end
 
 function count_sub_triangulation(m,mpoint_to_value)
-  n_scells = 0
-  n_spoints = 0
   mcell_to_mpoints = get_cell_to_points(m)
   table = get_lookup_table(m)
+  count_sub_triangulation(table,mcell_to_mpoints,mpoint_to_value)
+end
+
+function count_sub_triangulation(table,mcell_to_mpoints,mpoint_to_value)
+  n_scells = 0
+  n_spoints = 0
   for mcell in 1:length(mcell_to_mpoints)
     case = compute_case(mcell_to_mpoints,mpoint_to_value,mcell)
     n_scells += length(table.case_to_subcell_to_inout[case])
@@ -444,24 +448,18 @@ end
 
 function initial_sub_triangulation(grid::Grid,geom::DiscreteGeometry)
   ugrid = UnstructuredGrid(grid)
-  tree = get_tree(geom)
-  ls_to_point_to_value, oid_to_ls = _find_unique_leaves(tree)
+  ls_to_point_to_value, oid_to_ls = _find_unique_leaves(get_tree(geom))
   out = _initial_sub_triangulation(ugrid,ls_to_point_to_value)
   out[1], out[2], out[3], oid_to_ls
 end
 
 function _initial_sub_triangulation(grid::UnstructuredGrid,ls_to_point_to_value)
-
   cutgrid, ls_to_cutpoint_to_value, ls_to_bgcell_to_inoutcut = _extract_grid_of_cut_cells(grid,ls_to_point_to_value)
-
   subtrian, ls_to_subpoint_to_value = _simplexify_and_isolate_cells_in_cutgrid(cutgrid,ls_to_cutpoint_to_value)
-
-  subtrian, ls_to_subpoint_to_value, ls_to_bgcell_to_inoutcut
+  return subtrian, ls_to_subpoint_to_value, ls_to_bgcell_to_inoutcut
 end
 
 function _extract_grid_of_cut_cells(grid,ls_to_point_to_value)
-
-
   p = _check_and_get_polytope(grid)
   table = LookupTable(p)
   cell_to_points = get_cell_node_ids(grid)
@@ -682,4 +680,4 @@ function  _ensure_positive_jacobians_facets_work!(tcell_to_tpoints,c1,c2,tjac_q,
       tcell_to_tpoints.data[p+2] = p1
     end
   end
-end
\ No newline at end of file
+end

From c305cfb6099d875cd0d6d4396e56af2a17564437 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Wed, 5 Mar 2025 09:33:21 +1100
Subject: [PATCH 09/11] Update CI actions

---
 .github/dependabot.yml   |  7 ++++++
 .github/workflows/ci.yml | 49 +++++++++-------------------------------
 2 files changed, 18 insertions(+), 38 deletions(-)
 create mode 100644 .github/dependabot.yml

diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 00000000..d60f0707
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "monthly"
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index f316fd1c..441df67b 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -14,21 +14,12 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
@@ -48,21 +39,12 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - run: |
           julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
@@ -83,21 +65,12 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - run: |
           julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
@@ -109,8 +82,8 @@ jobs:
     name: Documentation
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1.8'
       - run: |

From 28d4404efe000407be018af9fbc36ab9b5186620 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Wed, 5 Mar 2025 09:37:52 +1100
Subject: [PATCH 10/11] Update x86 CI actions

---
 .github/workflows/ci_x86.yml | 15 +++------------
 1 file changed, 3 insertions(+), 12 deletions(-)

diff --git a/.github/workflows/ci_x86.yml b/.github/workflows/ci_x86.yml
index 03ca2ab1..4415f5a5 100644
--- a/.github/workflows/ci_x86.yml
+++ b/.github/workflows/ci_x86.yml
@@ -20,20 +20,11 @@ jobs:
         arch:
           - x86 
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1

From 921a96c156ca0cc6338aa366102b00f8e84abe91 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Wed, 5 Mar 2025 14:57:06 +1100
Subject: [PATCH 11/11] Added temp dirs

---
 .github/workflows/ci.yml                            |  8 ++++----
 test/AgFEMTests/PeriodicAgFEMSpacesTests.jl         |  3 ++-
 test/AlgoimUtilsTests/VisualizationTests.jl         |  5 +++--
 test/DistributedTests/AggregationTests.jl           |  9 ++++-----
 .../DistributedDiscreteGeometryPoissonTest.jl       |  7 ++++---
 .../DistributedLSDiscreteGeometryPoissonTest.jl     |  7 ++++---
 ...eriodicDistributedDiscreteGeometryPoissonTest.jl |  7 ++++---
 test/DistributedTests/PoissonTests.jl               |  7 ++++---
 test/DistributedTests/mpi/runtests.jl               | 10 +++++-----
 test/DistributedTests/testing_remote_no_aggs.jl     | 13 +++++++------
 .../EmbeddedBimaterialPoissonCutFEMTests.jl         |  5 +++--
 test/GridapEmbeddedTests/TraceFEMTests.jl           |  2 --
 12 files changed, 44 insertions(+), 39 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 441df67b..d116b290 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -8,7 +8,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.8'
+          - '1.10'
         os:
           - ubuntu-latest
         arch:
@@ -33,7 +33,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.8'
+          - '1.10'
         os:
           - ubuntu-latest
         arch:
@@ -59,7 +59,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.8'
+          - '1.10'
         os:
           - ubuntu-latest
         arch:
@@ -85,7 +85,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: '1.8'
+          version: '1.10'
       - run: |
           julia --project=docs -e '
             using Pkg
diff --git a/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl b/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl
index 26b34bae..9627e598 100644
--- a/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl
+++ b/test/AgFEMTests/PeriodicAgFEMSpacesTests.jl
@@ -40,7 +40,8 @@ U = TrialFESpace(Vagg)
 v(x) = (x[1]-0.5)^2 + (x[2]-0.5)^2
 vhagg = interpolate(v,Vagg)
 
-writevtk(Ω_ac,"test",cellfields=["v"=>vhagg])
+path = mktempdir()
+writevtk(Ω_ac,joinpath(path,"test"),cellfields=["v"=>vhagg])
 
 tol = 10e-7
 @test sum( ∫(abs2(v-vhagg))dΩ ) < tol
diff --git a/test/AlgoimUtilsTests/VisualizationTests.jl b/test/AlgoimUtilsTests/VisualizationTests.jl
index 91181789..d0b52930 100644
--- a/test/AlgoimUtilsTests/VisualizationTests.jl
+++ b/test/AlgoimUtilsTests/VisualizationTests.jl
@@ -37,7 +37,8 @@ module VisualizationTests
   vquad = Quadrature(algoim,phi,degree,phase=IN)
   _,dΩ = TriangulationAndMeasure(Ω,vquad)
 
-  writevtk(dΓ,"res_sur",cellfields=["f"=>fₕ],qhulltype=convexhull)
-  writevtk([dΩ,dΓ],"res_vol",cellfields=["f"=>fₕ])
+  path = mktempdir()
+  writevtk(dΓ,joinpath(path,"res_sur"),cellfields=["f"=>fₕ],qhulltype=convexhull)
+  writevtk([dΩ,dΓ],joinpath("res_vol"),cellfields=["f"=>fₕ])
 
 end # module
\ No newline at end of file
diff --git a/test/DistributedTests/AggregationTests.jl b/test/DistributedTests/AggregationTests.jl
index 0dff75cf..7175232b 100644
--- a/test/DistributedTests/AggregationTests.jl
+++ b/test/DistributedTests/AggregationTests.jl
@@ -38,8 +38,6 @@ bgf_to_ioc = compute_bgfacet_to_inoutcut(bgmodel,geo)
 
 Ω = Triangulation(cutgeo)
 
-writevtk(Ω,"trian")
-
 strategy = AggregateCutCellsByThreshold(1.0)
 aggregates,aggregate_owner,aggregate_neig = distributed_aggregate(
   strategy,cutgeo,geo,IN)
@@ -74,9 +72,10 @@ end
 Ωin = Triangulation(cutgeo,IN)
 Γ = EmbeddedBoundary(cutgeo)
 
-writevtk(Ωin,"trian_in")
-writevtk(Γ,"bnd")
-writevtk(Ωbg,"bgtrian",celldata=
+path = mktempdir()
+writevtk(Ωin,joinpath(path,"trian_in"))
+writevtk(Γ,joinpath(path,"bnd"))
+writevtk(Ωbg,joinpath(path,"bgtrian"),celldata=
    ["aggregate"=>oaggregates,
     "aggregate_owner"=>oaggregate_owner])
 
diff --git a/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl
index 10a99591..9388124d 100644
--- a/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl
+++ b/test/DistributedTests/DistributedDiscreteGeometryPoissonTest.jl
@@ -98,15 +98,16 @@ function main(distribute,parts;
     map(Reindex(col),oid)
   end
 
-  writevtk(Ω_bg,"trian",
+  path = mktempdir()
+  writevtk(Ω_bg,joinpath(path,"trian"),
     celldata=[
       "aggregate"=>own_aggregates,
       "color"=>own_colors,
       "gid"=>own_to_global(gids)])#,
   #  cellfields=["uh"=>uh])
 
-  writevtk(Ω,"trian_O",cellfields=["uh"=>uh])
-  writevtk(Γ,"trian_G")
+  writevtk(Ω,joinpath(path,"trian_O"),cellfields=["uh"=>uh])
+  writevtk(Γ,joinpath(path,"trian_G"))
   @test el2/ul2 < 1.e-8
   @test eh1/uh1 < 1.e-7
 
diff --git a/test/DistributedTests/DistributedLSDiscreteGeometryPoissonTest.jl b/test/DistributedTests/DistributedLSDiscreteGeometryPoissonTest.jl
index 68bbf39f..d27d8f7c 100644
--- a/test/DistributedTests/DistributedLSDiscreteGeometryPoissonTest.jl
+++ b/test/DistributedTests/DistributedLSDiscreteGeometryPoissonTest.jl
@@ -93,15 +93,16 @@ function main(distribute,parts;
     map(Reindex(col),oid)
   end
 
-  writevtk(Ω_bg,"trian",
+  path = mktempdir()
+  writevtk(Ω_bg,joinpath(path,"trian"),
     celldata=[
       "aggregate"=>own_aggregates,
       "color"=>own_colors,
       "gid"=>own_to_global(gids)])#,
   #  cellfields=["uh"=>uh])
 
-  writevtk(Ω,"trian_O",cellfields=["uh"=>uh])
-  writevtk(Γ,"trian_G")
+  writevtk(Ω,joinpath(path,"trian_O"),cellfields=["uh"=>uh])
+  writevtk(Γ,joinpath(path,"trian_G"))
   @test el2/ul2 < 1.e-8
   @test eh1/uh1 < 1.e-7
 
diff --git a/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl b/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl
index 26cde8c3..2c7de035 100644
--- a/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl
+++ b/test/DistributedTests/PeriodicDistributedDiscreteGeometryPoissonTest.jl
@@ -151,15 +151,16 @@ function main(distribute,parts;
     map(Reindex(col),oid)
   end
 
-  writevtk(Ω_bg,"trian",
+  path = mktempdir()
+  writevtk(Ω_bg,joinpath(path,"trian"),
     celldata=[
       "aggregate"=>own_aggregates,
       "color"=>own_colors,
       "gid"=>own_to_global(gids)])#,
   #  cellfields=["uh"=>uh])
 
-  writevtk(Ω,"trian_O",cellfields=["uh"=>uh])
-  writevtk(Γ,"trian_G")
+  writevtk(Ω,joinpath(path,"trian_O"),cellfields=["uh"=>uh])
+  writevtk(Γ,joinpath(path,"trian_G"))
   @test el2/ul2 < 1.e-8
   @test eh1/uh1 < 1.e-7
 
diff --git a/test/DistributedTests/PoissonTests.jl b/test/DistributedTests/PoissonTests.jl
index feb6e9b2..d5ece5b9 100644
--- a/test/DistributedTests/PoissonTests.jl
+++ b/test/DistributedTests/PoissonTests.jl
@@ -96,15 +96,16 @@ function main(distribute,parts;
     map(Reindex(col),oid)
   end
 
-  writevtk(Ω_bg,"trian",
+  path = mktempdir()
+  writevtk(Ω_bg,joinpath(path,"trian"),
     celldata=[
       "aggregate"=>own_aggregates,
       "color"=>own_colors,
       "gid"=>own_to_global(gids)])#,
   #  cellfields=["uh"=>uh])
 
-  writevtk(Ω,"trian_O",cellfields=["uh"=>uh,"eh"=>e])
-  writevtk(Γ,"trian_G")
+  writevtk(Ω,joinpath(path,"trian_O"),cellfields=["uh"=>uh,"eh"=>e])
+  writevtk(Γ,joinpath(path,"trian_G"))
   @test el2/ul2 < 1.e-8
   @test eh1/uh1 < 1.e-7
 
diff --git a/test/DistributedTests/mpi/runtests.jl b/test/DistributedTests/mpi/runtests.jl
index 4f2c40c9..44a1968b 100644
--- a/test/DistributedTests/mpi/runtests.jl
+++ b/test/DistributedTests/mpi/runtests.jl
@@ -6,8 +6,8 @@ using MPI
 #Sysimage
 sysimage=nothing
 if length(ARGS)==1
-   @assert isfile(ARGS[1]) "$(ARGS[1]) must be a valid Julia sysimage file"
-   sysimage=ARGS[1]
+  @assert isfile(ARGS[1]) "$(ARGS[1]) must be a valid Julia sysimage file"
+  sysimage=ARGS[1]
 end
 
 mpidir = @__DIR__
@@ -17,10 +17,10 @@ repodir = joinpath(testdir,"..","..")
 function run_driver(procs,file,sysimage)
   mpiexec() do cmd
     if sysimage!=nothing
-       extra_args="-J$(sysimage)"
-       run(`$cmd -n $procs $(Base.julia_cmd()) $(extra_args) --project=$repodir $(joinpath(mpidir,file))`)
+      extra_args="-J$(sysimage)"
+      run(`$cmd -n $procs $(Base.julia_cmd()) $(extra_args) --project=$repodir $(joinpath(mpidir,file))`)
     else
-       run(`$cmd -n $procs $(Base.julia_cmd()) --project=$repodir $(joinpath(mpidir,file))`)
+      run(`$cmd -n $procs $(Base.julia_cmd()) --project=$repodir $(joinpath(mpidir,file))`)
     end
     @test true
   end
diff --git a/test/DistributedTests/testing_remote_no_aggs.jl b/test/DistributedTests/testing_remote_no_aggs.jl
index 87e5d1aa..0bee27bd 100644
--- a/test/DistributedTests/testing_remote_no_aggs.jl
+++ b/test/DistributedTests/testing_remote_no_aggs.jl
@@ -133,18 +133,19 @@ uh1 = h1(uh)
 Γ = EmbeddedBoundary(cutgeo)
 Ωbg = Triangulation(bgmodel)
 
-writevtk(Ω,"trian");
-writevtk(Γ,"bnd");
-writevtk(Ωbg,"bg_trian");
-writevtk(Ω_act,"act_trian");
+path = mktempdir()
+writevtk(Ω,joinpath(path,"trian"));
+writevtk(Γ,joinpath(path,"bnd"));
+writevtk(Ωbg,joinpath(path,"bg_trian"));
+writevtk(Ω_act,joinpath(path,"act_trian"));
 
-writevtk(Ω,"trian",
+writevtk(Ω,joinpath(path,"trian"),
   cellfields=["uh"=>uh,"u"=>u,"e"=>e],);
 
 
 map(local_views(uh),local_views(bgmodel),ranks) do uh,m,p
   trian = Triangulation(m)
-  writevtk(trian,"ltrian_$p",cellfields=["uh"=>uh])
+  writevtk(trian,joinpath(path,"ltrian_$p"),cellfields=["uh"=>uh])
 end
 
 end # module
diff --git a/test/GridapEmbeddedTests/EmbeddedBimaterialPoissonCutFEMTests.jl b/test/GridapEmbeddedTests/EmbeddedBimaterialPoissonCutFEMTests.jl
index 3d3e718d..138857b5 100644
--- a/test/GridapEmbeddedTests/EmbeddedBimaterialPoissonCutFEMTests.jl
+++ b/test/GridapEmbeddedTests/EmbeddedBimaterialPoissonCutFEMTests.jl
@@ -103,10 +103,11 @@ uh1, uh2 = solve(op)
 uh = (uh1,uh2)
 
 # Postprocess
+path = mktempdir()
 qh1 = α1*∇(uh1)
 qh2 = α2*∇(uh2)
-writevtk(Ω1,"results1",cellfields=["uh"=>uh1,"qh"=>qh1])
-writevtk(Ω2,"results2",cellfields=["uh"=>uh2,"qh"=>qh2])
+writevtk(Ω1,joinpath(path,"results1"),cellfields=["uh"=>uh1,"qh"=>qh1])
+writevtk(Ω2,joinpath(path,"results2"),cellfields=["uh"=>uh2,"qh"=>qh2])
 
 #writevtk(model1,"model1")
 #writevtk(model2,"model2")
diff --git a/test/GridapEmbeddedTests/TraceFEMTests.jl b/test/GridapEmbeddedTests/TraceFEMTests.jl
index 7953e8f0..0e46c7d2 100644
--- a/test/GridapEmbeddedTests/TraceFEMTests.jl
+++ b/test/GridapEmbeddedTests/TraceFEMTests.jl
@@ -26,8 +26,6 @@ cutgeom = cut(bgmodel,geom)
 Γ = EmbeddedBoundary(cutgeom,geom)
 Γg = GhostSkeleton(cutgeom,CUT,geom)
 
-writevtk(Γg,"Γg")
-
 order=1
 V = TestFESpace(Ωc,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
 U = TrialFESpace(V)