From 610d3621a14113c6f202a48737f6d168b6415fca Mon Sep 17 00:00:00 2001
From: zjwegert <60646897+zjwegert@users.noreply.github.com>
Date: Fri, 24 Oct 2025 16:45:27 +1000
Subject: [PATCH 1/7] Fix GridapEmbedded#115

---
 src/Interfaces/EmbeddedDiscretizations.jl | 40 +++++++++----
 test/IssuesTests/GridapEmbedded#115.jl    | 68 +++++++++++++++++++++++
 test/IssuesTests/runtests.jl              |  7 +++
 test/runtests.jl                          |  2 +
 4 files changed, 107 insertions(+), 10 deletions(-)
 create mode 100644 test/IssuesTests/GridapEmbedded#115.jl
 create mode 100644 test/IssuesTests/runtests.jl

diff --git a/src/Interfaces/EmbeddedDiscretizations.jl b/src/Interfaces/EmbeddedDiscretizations.jl
index 82c9486b..62615fd1 100644
--- a/src/Interfaces/EmbeddedDiscretizations.jl
+++ b/src/Interfaces/EmbeddedDiscretizations.jl
@@ -6,7 +6,7 @@ abstract type AbstractEmbeddedDiscretization <: GridapType end
 """
     struct EmbeddedDiscretization{Dc,T} <: AbstractEmbeddedDiscretization
 
-This structure contains all the required information to build integration `Triangulation`s 
+This structure contains all the required information to build integration `Triangulation`s
 for a cut model.
 
 ## Constructors
@@ -19,11 +19,11 @@ for a cut model.
 - `geo::CSG.Geometry`: the geometry used to cut the background mesh
 - `subcells::SubCellData`: collection of cut subcells, attached to the background mesh
 - `subfacets::SubFacetData`: collection of cut facets, attached to the background mesh
-- `ls_to_bgcell_to_inoutcut::Vector{Vector{Int8}}`: list of IN/OUT/CUT states for each cell 
+- `ls_to_bgcell_to_inoutcut::Vector{Vector{Int8}}`: list of IN/OUT/CUT states for each cell
    in the background mesh, for each node in the geometry tree.
-- `ls_to_subcell_to_inoutcut::Vector{Vector{Int8}}`: list of IN/OUT/CUT states for each subcell 
+- `ls_to_subcell_to_inoutcut::Vector{Vector{Int8}}`: list of IN/OUT/CUT states for each subcell
    in the cut part of the mesh, for each node in the geometry tree.
-- `ls_to_subfacet_to_inoutcut::Vector{Vector{Int8}}`: list of IN/OUT/CUT states for each subfacet 
+- `ls_to_subfacet_to_inoutcut::Vector{Vector{Int8}}`: list of IN/OUT/CUT states for each subfacet
    in the cut part of the mesh, for each node in the geometry tree.
 
 ## Methods
@@ -259,12 +259,12 @@ end
 """
     Triangulation(cut::EmbeddedDiscretization[,in_or_out=PHYSICAL_IN])
 
-Creates a triangulation containing the cell and subcells of the embedded domain selected by 
+Creates a triangulation containing the cell and subcells of the embedded domain selected by
 `in_or_out`.
 
 - If only background cells are selected, the result will be a regular Gridap triangulation.
 - If only subcells are selected, the result will be a [`SubCellTriangulation`](@ref).
-- If both background cells and subcells are selected, the result will be an `AppendedTriangulation`, 
+- If both background cells and subcells are selected, the result will be an `AppendedTriangulation`,
   containing a [`SubCellTriangulation`](@ref) and a regular Gridap triangulation.
 
 """
@@ -317,11 +317,21 @@ function Triangulation(cut::EmbeddedDiscretization,in_or_out::Integer,geo::CSG.G
 end
 
 function Triangulation(cut::EmbeddedDiscretization,in_or_out::ActiveInOrOut,geo::CSG.Geometry)
-  cell_to_inoutcut = compute_bgcell_to_inoutcut(cut,geo)
-  cell_to_mask = lazy_map(i-> i==CUT || i==in_or_out.in_or_out,cell_to_inoutcut)
+  cell_to_mask = compute_mask(cut,in_or_out,geo)
   Triangulation(cut.bgmodel,cell_to_mask)
 end
 
+function compute_mask(cut::EmbeddedDiscretization,in_or_out::Union{CutInOrOut,ActiveInOrOut},geo::CSG.Geometry)
+  bgcell_to_inoutcut = compute_bgcell_to_inoutcut(cut,geo)
+  subcell_to_inoutcut = lazy_map(Reindex(bgcell_to_inoutcut),cut.subcells.cell_to_bgcell)
+  subcell_to_inout = compute_subcell_to_inout(cut,geo)
+  mask = lazy_map( (a,b) -> a==CUT && b==in_or_out.in_or_out, subcell_to_inoutcut, subcell_to_inout   )
+  newsubcells = findall(mask)
+  cell_to_mask = collect(Bool,bgcell_to_inoutcut .== in_or_out.in_or_out) # in/out cells
+  cell_to_mask[cut.subcells.cell_to_bgcell[newsubcells]] .= true # (true) cut cells
+  return cell_to_mask
+end
+
 function compute_subcell_to_inout(cut::EmbeddedDiscretization,geo::CSG.Geometry)
 
   tree = get_tree(geo)
@@ -465,10 +475,10 @@ end
 """
     GhostSkeleton(cut::EmbeddedDiscretization[,in_or_out=ACTIVE_IN])
 
-Creates a triangulation containing the ghost facets. Ghosts facets are defined as the facets 
+Creates a triangulation containing the ghost facets. Ghosts facets are defined as the facets
 of the **background mesh** that are adjacent to at least one `CUT` background cell.
 
-Mostly used for CUT-FEM stabilisation. 
+Mostly used for CUT-FEM stabilisation.
 """
 function GhostSkeleton(cut::EmbeddedDiscretization)
   GhostSkeleton(cut,ACTIVE_IN)
@@ -506,6 +516,16 @@ function GhostSkeleton(cut::EmbeddedDiscretization,in_or_out,geo::CSG.Geometry)
   @notimplementedif in_or_out in (PHYSICAL_IN,PHYSICAL_OUT) "Not implemented but not needed in practice. Ghost stabilization can be integrated in full facets."
   @assert in_or_out in (ACTIVE_IN,ACTIVE_OUT) || in_or_out.in_or_out in (IN,OUT,CUT)
   cell_to_inoutcut = compute_bgcell_to_inoutcut(cut,geo)
+  cell_to_mask = compute_mask(cut,in_or_out,geo)
+  cut_bgcells = findall(x->x==CUT,cell_to_inoutcut)
+  # Check if cell is not truely cut
+  for c in cut_bgcells
+      if !cell_to_mask[c]
+        # It doesn't matter if they're marked as in or out for purpose of GhostSkeleton
+        cell_to_inoutcut[c] = IN
+      end
+  end
+  # build facet mask
   model = cut.bgmodel
   topo = get_grid_topology(model)
   D = num_cell_dims(model)
diff --git a/test/IssuesTests/GridapEmbedded#115.jl b/test/IssuesTests/GridapEmbedded#115.jl
new file mode 100644
index 00000000..3705a1e4
--- /dev/null
+++ b/test/IssuesTests/GridapEmbedded#115.jl
@@ -0,0 +1,68 @@
+using Gridap, GridapEmbedded
+using Gridap.Geometry, Gridap.Adaptivity
+using GridapEmbedded.LevelSetCutters
+using GridapEmbedded.LevelSetCutters: SubCellTriangulation
+using Test
+
+### Test 1
+n = 12
+base_model = UnstructuredDiscreteModel(CartesianDiscreteModel((0,1,0,1),(n,n)))
+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))
+
+Ω2 = Triangulation(cutgeo,PHYSICAL,"Ω2")
+Ω2_act = Triangulation(cutgeo,ACTIVE,"Ω2")
+Λ_Ω2 = GhostSkeleton(cutgeo,"Ω2")
+
+@test iszero(num_cells(Ω2))
+@test iszero(num_cells(Ω2_act))
+@test iszero(num_cells(Λ_Ω2))
+
+### Test 2
+n = 48
+base_model = UnstructuredDiscreteModel(CartesianDiscreteModel((0,1,0,1),(n,n)))
+ref_model = refine(base_model, refinement_method = "barycentric")
+model = get_model(ref_model)
+f((x,y),a1,a2,b1,b2,c) = cos(a1*pi*(x-a2))*cos(b1*pi*(y-b2)) - c
+φ1((x,y)) = f((x,y),4,1/8,4,1/8,0.5)
+φ2((x,y)) = f((x,y),4,3/8,4,1/8,0.5)
+φ3((x,y)) = f((x,y),4,3/8,4,1/8,0.5-0.5/n)
+V_φ = TestFESpace(model,ReferenceFE(lagrangian,Float64,1))
+φh1 = interpolate(φ1,V_φ); φh2 = interpolate(φ2,V_φ); φh3 = interpolate(φ3,V_φ)
+
+geo1 = DiscreteGeometry(φh1,model,name="φ1")
+geo2 = DiscreteGeometry(φh2,model,name="φ2")
+geo3 = DiscreteGeometry(φh3,model,name="φ3")
+Ω1_geo = intersect(geo1,geo2,name="Ω1")
+Ω2_geo = setdiff(geo3,Ω1_geo,name="Ω2")
+cutgeo = cut(model,union(Ω1_geo,Ω2_geo))
+
+Ω2 = Triangulation(cutgeo,PHYSICAL,"Ω2")
+Ω2_act = Triangulation(cutgeo,ACTIVE,"Ω2")
+Λ_Ω2 = GhostSkeleton(cutgeo,"Ω2")
+
+# Check cells are in subcell to bgcell list or in cell to bgcell list
+for c in Ω2_act.tface_to_mface
+  @test c ∈ Ω2.a.subcells.cell_to_bgcell || c ∈ Ω2.b.tface_to_mface
+end
+for c in Λ_Ω2.plus.glue.face_to_cell
+  @test c ∈ Ω2.a.subcells.cell_to_bgcell || c ∈ Ω2.b.tface_to_mface
+end
+for c in Λ_Ω2.minus.glue.face_to_cell
+  @test c ∈ Ω2.a.subcells.cell_to_bgcell || c ∈ Ω2.b.tface_to_mface
+end
+
+# writevtk(Triangulation(model),"results/Omega")
+# writevtk(Triangulation(cutgeo,"Ω2"),"results/Omega2")
+# writevtk(Triangulation(cutgeo,ACTIVE,"Ω2"),"results/Omega2_act")
+# writevtk(Λ_Ω2,"results/Omega2_ghost_skel")
\ No newline at end of file
diff --git a/test/IssuesTests/runtests.jl b/test/IssuesTests/runtests.jl
new file mode 100644
index 00000000..56088007
--- /dev/null
+++ b/test/IssuesTests/runtests.jl
@@ -0,0 +1,7 @@
+module IssuesTests
+
+using Test
+
+@testset "GridapEmbedded#115" begin include("GridapEmbedded#115.jl") end
+
+end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index 31fc7176..7651185f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,6 +23,8 @@ if QHULL_LOADED
 
   @time @testset "Distributed" begin include("DistributedTests/runtests.jl") end
 
+  @time @testset "Issues" begin include("IssuesTests/runtests.jl") end
+
   include(joinpath(@__DIR__,"..","examples","runexamples.jl"))
 
 else

From 228c37eb200f5edc13c7a2955eb15bdd9e1e8b55 Mon Sep 17 00:00:00 2001
From: zjwegert <60646897+zjwegert@users.noreply.github.com>
Date: Fri, 24 Oct 2025 16:46:44 +1000
Subject: [PATCH 2/7] Update test

---
 test/IssuesTests/GridapEmbedded#115.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/test/IssuesTests/GridapEmbedded#115.jl b/test/IssuesTests/GridapEmbedded#115.jl
index 3705a1e4..21457e06 100644
--- a/test/IssuesTests/GridapEmbedded#115.jl
+++ b/test/IssuesTests/GridapEmbedded#115.jl
@@ -1,3 +1,4 @@
+module GridapEmbeddedIssue115Tests
 using Gridap, GridapEmbedded
 using Gridap.Geometry, Gridap.Adaptivity
 using GridapEmbedded.LevelSetCutters
@@ -65,4 +66,5 @@ end
 # writevtk(Triangulation(model),"results/Omega")
 # writevtk(Triangulation(cutgeo,"Ω2"),"results/Omega2")
 # writevtk(Triangulation(cutgeo,ACTIVE,"Ω2"),"results/Omega2_act")
-# writevtk(Λ_Ω2,"results/Omega2_ghost_skel")
\ No newline at end of file
+# writevtk(Λ_Ω2,"results/Omega2_ghost_skel")
+end
\ No newline at end of file

From 1e528e774977c611e967853a6f947e2c4450e11a Mon Sep 17 00:00:00 2001
From: zjwegert <60646897+zjwegert@users.noreply.github.com>
Date: Fri, 24 Oct 2025 17:12:31 +1000
Subject: [PATCH 3/7] Conform to Gridap tests

---
 .../issues/issue_115.jl}                               | 10 ++++------
 test/GridapEmbeddedTests/runtests.jl                   |  4 ++++
 test/IssuesTests/runtests.jl                           |  7 -------
 test/runtests.jl                                       |  2 --
 4 files changed, 8 insertions(+), 15 deletions(-)
 rename test/{IssuesTests/GridapEmbedded#115.jl => GridapEmbeddedTests/issues/issue_115.jl} (87%)
 delete mode 100644 test/IssuesTests/runtests.jl

diff --git a/test/IssuesTests/GridapEmbedded#115.jl b/test/GridapEmbeddedTests/issues/issue_115.jl
similarity index 87%
rename from test/IssuesTests/GridapEmbedded#115.jl
rename to test/GridapEmbeddedTests/issues/issue_115.jl
index 21457e06..13a15774 100644
--- a/test/IssuesTests/GridapEmbedded#115.jl
+++ b/test/GridapEmbeddedTests/issues/issue_115.jl
@@ -1,4 +1,3 @@
-module GridapEmbeddedIssue115Tests
 using Gridap, GridapEmbedded
 using Gridap.Geometry, Gridap.Adaptivity
 using GridapEmbedded.LevelSetCutters
@@ -63,8 +62,7 @@ for c in Λ_Ω2.minus.glue.face_to_cell
   @test c ∈ Ω2.a.subcells.cell_to_bgcell || c ∈ Ω2.b.tface_to_mface
 end
 
-# writevtk(Triangulation(model),"results/Omega")
-# writevtk(Triangulation(cutgeo,"Ω2"),"results/Omega2")
-# writevtk(Triangulation(cutgeo,ACTIVE,"Ω2"),"results/Omega2_act")
-# writevtk(Λ_Ω2,"results/Omega2_ghost_skel")
-end
\ No newline at end of file
+writevtk(Triangulation(model),"Omega")
+writevtk(Triangulation(cutgeo,"Ω2"),"Omega2")
+writevtk(Triangulation(cutgeo,ACTIVE,"Ω2"),"Omega2_act")
+writevtk(Λ_Ω2,"Omega2_ghost_skel")
\ No newline at end of file
diff --git a/test/GridapEmbeddedTests/runtests.jl b/test/GridapEmbeddedTests/runtests.jl
index f6d52df3..2c89ead7 100644
--- a/test/GridapEmbeddedTests/runtests.jl
+++ b/test/GridapEmbeddedTests/runtests.jl
@@ -20,4 +20,8 @@ using Test
 
 @time @testset "TraceFEM" begin include("TraceFEMTests.jl") end
 
+@testset "Issue reproducers" begin
+  @time @testset "Issue614"  begin include("issues/issue_115.jl")  end
+end
+
 end
diff --git a/test/IssuesTests/runtests.jl b/test/IssuesTests/runtests.jl
deleted file mode 100644
index 56088007..00000000
--- a/test/IssuesTests/runtests.jl
+++ /dev/null
@@ -1,7 +0,0 @@
-module IssuesTests
-
-using Test
-
-@testset "GridapEmbedded#115" begin include("GridapEmbedded#115.jl") end
-
-end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index 7651185f..31fc7176 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,8 +23,6 @@ if QHULL_LOADED
 
   @time @testset "Distributed" begin include("DistributedTests/runtests.jl") end
 
-  @time @testset "Issues" begin include("IssuesTests/runtests.jl") end
-
   include(joinpath(@__DIR__,"..","examples","runexamples.jl"))
 
 else

From e6482172bd565ef6724877d4f2c04d2b854c5057 Mon Sep 17 00:00:00 2001
From: zjwegert <60646897+zjwegert@users.noreply.github.com>
Date: Fri, 24 Oct 2025 17:19:05 +1000
Subject: [PATCH 4/7] Fix dependency

---
 test/IssuesTests/runtests.jl | 7 +++++++
 1 file changed, 7 insertions(+)
 create mode 100644 test/IssuesTests/runtests.jl

diff --git a/test/IssuesTests/runtests.jl b/test/IssuesTests/runtests.jl
new file mode 100644
index 00000000..56088007
--- /dev/null
+++ b/test/IssuesTests/runtests.jl
@@ -0,0 +1,7 @@
+module IssuesTests
+
+using Test
+
+@testset "GridapEmbedded#115" begin include("GridapEmbedded#115.jl") end
+
+end # module

From 4c3049067c115b549d3ce422d516fa231ba08625 Mon Sep 17 00:00:00 2001
From: zjwegert <60646897+zjwegert@users.noreply.github.com>
Date: Fri, 24 Oct 2025 17:19:54 +1000
Subject: [PATCH 5/7] Undo, fix dep

---
 Project.toml                 | 2 +-
 test/IssuesTests/runtests.jl | 7 -------
 2 files changed, 1 insertion(+), 8 deletions(-)
 delete mode 100644 test/IssuesTests/runtests.jl

diff --git a/Project.toml b/Project.toml
index b5e91a35..5c8452a7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -29,7 +29,7 @@ CxxWrap = "0.16"
 FillArrays = "0.10, 0.11, 0.12, 0.13, 1"
 FiniteDiff = "2.27.0"
 ForwardDiff = "0.10.38, 1"
-Graphs = "1.12.0"
+Graphs = "1.13"
 Gridap = "0.18.12, 0.19"
 GridapDistributed = "0.3, 0.4"
 MPI = "0.20"
diff --git a/test/IssuesTests/runtests.jl b/test/IssuesTests/runtests.jl
deleted file mode 100644
index 56088007..00000000
--- a/test/IssuesTests/runtests.jl
+++ /dev/null
@@ -1,7 +0,0 @@
-module IssuesTests
-
-using Test
-
-@testset "GridapEmbedded#115" begin include("GridapEmbedded#115.jl") end
-
-end # module

From cec20ff65a137c273bbae5bcf742d1472b980e38 Mon Sep 17 00:00:00 2001
From: zjwegert <60646897+zjwegert@users.noreply.github.com>
Date: Fri, 24 Oct 2025 17:35:36 +1000
Subject: [PATCH 6/7] Fix deps

---
 Project.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/Project.toml b/Project.toml
index 5c8452a7..d4a47936 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 Algoim = "0eb9048c-21de-4c7a-bfac-056de1940b74"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
@@ -26,6 +27,7 @@ AbstractTrees = "0.3.3, 0.4"
 Algoim = "0.2.2"
 Combinatorics = "1"
 CxxWrap = "0.16"
+DataStructures = "0.18"
 FillArrays = "0.10, 0.11, 0.12, 0.13, 1"
 FiniteDiff = "2.27.0"
 ForwardDiff = "0.10.38, 1"

From 18a820287e3409a651bac3829633da9a2a2a01d5 Mon Sep 17 00:00:00 2001
From: Zachary J Wegert <60646897+zjwegert@users.noreply.github.com>
Date: Thu, 13 Nov 2025 08:45:06 +1000
Subject: [PATCH 7/7] Update DataStructures version to 0.19

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index d4a47936..014bb2fe 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,7 +27,7 @@ AbstractTrees = "0.3.3, 0.4"
 Algoim = "0.2.2"
 Combinatorics = "1"
 CxxWrap = "0.16"
-DataStructures = "0.18"
+DataStructures = "0.19"
 FillArrays = "0.10, 0.11, 0.12, 0.13, 1"
 FiniteDiff = "2.27.0"
 ForwardDiff = "0.10.38, 1"