Active triangulation can produce isolated elements

[zjwegert]

Hi all,
For context, this is an edge case that I've encountered with moving geometries in which one of the phases vanish.

The issue presents in the following setup: suppose we have two level-set functions that we use to define two distinct domains Ω1 and Ω2. For example, suppose Ω1 and Ω2 are the red and green domains in this image:

In grey, we visualise the ACTIVE domain for Ω2. This is produced by the following:

using Gridap, GridapEmbedded
using Gridap.Geometry, Gridap.Adaptivity
using GridapEmbedded.LevelSetCutters
using GridapEmbedded.LevelSetCutters: SubCellTriangulation
n = 12;
base_model = UnstructuredDiscreteModel(CartesianDiscreteModel((0,1,0,1),(12,12)))
ref_model = refine(base_model, refinement_method = "barycentric")
model = get_model(ref_model)
φ1((x,y)) = (x-0.5)^2 + (y-0.5)^2 - 0.3^2
φ2((x,y)) = (x-0.5)^2 + (y-0.5)^2 - (0.3+0.5/n)^2
V_φ = TestFESpace(model,ReferenceFE(lagrangian,Float64,1))
φh1 = interpolate(φ1,V_φ); φh2 = interpolate(φ2,V_φ)

geo1 = DiscreteGeometry(φh1,model,name="φ1")
geo2 = DiscreteGeometry(φh2,model,name="φ2")
Ω1_geo = setdiff(!geo1,geo2,name="Ω1")
Ω2_geo = setdiff(geo2,geo1,name="Ω2")
cutgeo = cut(model,union(Ω1_geo,Ω2_geo))

writevtk(Triangulation(cutgeo,"Ω1"),"results/Omega1")
writevtk(Triangulation(cutgeo,"Ω2"),"results/Omega2")
writevtk(Triangulation(cutgeo,ACTIVE,"Ω1"),"results/Omega1_act")
writevtk(Triangulation(cutgeo,ACTIVE,"Ω2"),"results/Omega2_act")
writevtk(EmbeddedBoundary(cutgeo,"Ω1","Ω2"),"results/Gamma12")

The physical geometry defining Ω2 vanishes if we change φ2 to φ2((x,y)) = (x-0.5)^2 + (y-0.5)^2 - (0.3-0.5/n)^2. However, the ACTIVE triangulation for Ω2 still has elements:

Depending on the variational problem we're solving, this can lead to elements that are isolated from the rest of the ACTIVE domain, resulting in a singular problem.

The above appears to happen because we compute the mask for the ACTIVE triangulation based on the maps from the background cells to IN/OUT/CUT and in the case of taking a setdiff, we mark overlapping CUT cells as CUT. This is almost always fine, except for the case above...

To get around this issue, I've been defining the ACTIVE triangulation based on the PHYSICAL triangulation (all background cells that are touched). I'm not sure this is an appropriate fix at a package level though...

function ActiveTriangulation(model::DiscreteModel,ttrian)
  mask = get_active_mask(model,ttrian)
  return Triangulation(model,mask)
end

function get_active_mask(model,t::AppendedTriangulation)
  mask = falses(num_cells(model))
  get_active_mask!(mask,t.a)
  get_active_mask!(mask,t.b)
  return mask
end

function get_active_mask(model,t::BodyFittedTriangulation)
  mask = falses(num_cells(model))
  get_active_mask!(mask,t)
  return mask
end

function get_active_mask!(mask::BitVector,t::SubCellTriangulation)
  mask[t.subcells.cell_to_bgcell] .= true
  nothing
end

function get_active_mask!(mask::BitVector,t::BodyFittedTriangulation)
  mask[t.tface_to_mface] .= true
  nothing
end

Edit: This also happens for GhostSkeleton