Merge branch 'master' of github.com:gridap/GridapEmbedded.jl into hdiv-stabilisation

Author: JordiManyer

.github/workflows/ci.yml

@@ -15,15 +15,15 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v6
-      - uses: julia-actions/setup-julia@v2
+      - uses: julia-actions/setup-julia@v3
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v2
+      - uses: julia-actions/cache@v3
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v5
+      - uses: codecov/codecov-action@v6
         with:
           file: lcov.info
   drivers:
@@ -40,11 +40,11 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v6
-      - uses: julia-actions/setup-julia@v2
+      - uses: julia-actions/setup-julia@v3
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v2
+      - uses: julia-actions/cache@v3
       - uses: julia-actions/julia-buildpkg@v1
       - run: |
           julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
@@ -66,11 +66,11 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v6
-      - uses: julia-actions/setup-julia@v2
+      - uses: julia-actions/setup-julia@v3
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v2
+      - uses: julia-actions/cache@v3
       - uses: julia-actions/julia-buildpkg@v1
       - run: |
           julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
@@ -83,7 +83,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v6
-      - uses: julia-actions/setup-julia@v2
+      - uses: julia-actions/setup-julia@v3
         with:
           version: '1.10'
       - run: |

.github/workflows/ci_x86.yml

@@ -21,10 +21,10 @@ jobs:
           - x86 
     steps:
       - uses: actions/checkout@v6
-      - uses: julia-actions/setup-julia@v2
+      - uses: julia-actions/setup-julia@v3
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v2
+      - uses: julia-actions/cache@v3
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1

NEWS.md

@@ -6,9 +6,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Fixed
+
+- Fix for issue #115, which is a bug in the logic of the `ACTIVE` mask and `GhostSkeleton` construction. Since PR [#123](https://github.com/gridap/GridapEmbedded.jl/pull/123).
+
 ## [0.9.11] - 2026-4-28
 
-- Small bugfix for empty embedded facet triangulations, following changes in Gridap v0.20.4. Since PR[#121](https://github.com/gridap/GridapEmbedded.jl/pull/121). 
+### Fixed
+
+- Small bugfix for empty embedded facet triangulations, following changes in Gridap v0.20.4. Since PR[#121](https://github.com/gridap/GridapEmbedded.jl/pull/121).
 
 ## [0.9.10] - 2026-3-24
 

src/Interfaces/EmbeddedDiscretizations.jl

@@ -316,9 +316,42 @@ function Triangulation(cut::EmbeddedDiscretization,in_or_out::Integer,geo::CSG.G
   Triangulation(cut.bgmodel,cell_to_mask)
 end
 
+function compute_bgcell_to_active_inoutcut(
+  cut::EmbeddedDiscretization, in_or_out, geo::CSG.Geometry
+)
+  compute_bgcell_to_active_inoutcut(cut, in_or_out.in_or_out, geo)
+end
+
+# Deals with the phantom cut issue (see Issue #115)
+function compute_bgcell_to_active_inoutcut(
+  cut::EmbeddedDiscretization, in_or_out::Integer, geo::CSG.Geometry
+)
+  @check in_or_out in (IN,OUT,CUT)
+  
+  bgcell_to_inoutcut = compute_bgcell_to_inoutcut(cut, geo)
+  if in_or_out == CUT
+    # CUT cells are always genuinely cut (no phantom issue), so no filtering needed
+    return bgcell_to_inoutcut
+  end
+
+  subcell_to_inout  = compute_subcell_to_inout(cut, geo)
+  subcell_to_bgcell = cut.subcells.cell_to_bgcell
+
+  keep = isequal(in_or_out)
+  bgcell_to_true_cut = fill(false, length(bgcell_to_inoutcut))
+  for (bgcell,io) in zip(subcell_to_bgcell, subcell_to_inout)
+    bgcell_to_true_cut[bgcell] |= keep(io)
+  end
+
+  out = -Int8(in_or_out) # IN → OUT, OUT → IN
+  map(bgcell_to_inoutcut, bgcell_to_true_cut) do ioc, is_true_cut
+    ifelse(ioc == CUT, ifelse(is_true_cut, ioc, out), ioc)
+  end
+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_inoutcut = compute_bgcell_to_active_inoutcut(cut, in_or_out, geo)
+  cell_to_mask = map(i -> i==CUT || i==in_or_out.in_or_out, cell_to_inoutcut)
   Triangulation(cut.bgmodel,cell_to_mask)
 end
 
@@ -486,40 +519,37 @@ function GhostSkeleton(cut::EmbeddedDiscretization,geo::CSG.Geometry)
   GhostSkeleton(cut,ACTIVE_IN,geo)
 end
 
-#function GhostSkeleton(cut::EmbeddedDiscretization,name::String,in_or_out)
-#  GhostSkeleton(cut,in_or_out,name)
-#end
-#
-#function GhostSkeleton(cut::EmbeddedDiscretization,geo::CSG.Geometry,in_or_out)
-#  GhostSkeleton(cut,in_or_out,geo)
-#end
-
 function GhostSkeleton(cut::EmbeddedDiscretization,in_or_out,name::String)
   geo = get_geometry(cut.geo,name)
   GhostSkeleton(cut,in_or_out,geo)
 end
 
-GhostSkeleton(cut::EmbeddedDiscretization,in_or_out::Integer,geo::CSG.Geometry) = GhostSkeleton(cut,ActiveInOrOut(in_or_out),geo)
+function GhostSkeleton(cut::EmbeddedDiscretization,in_or_out::Integer,geo::CSG.Geometry) 
+  GhostSkeleton(cut,ActiveInOrOut(in_or_out),geo)
+end
 
 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."
+  @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)
+
   model = cut.bgmodel
   topo = get_grid_topology(model)
   D = num_cell_dims(model)
+
+  cell_to_inoutcut = compute_bgcell_to_active_inoutcut(cut, in_or_out, geo)
   facet_to_cells = Table(get_faces(topo,D-1,D))
   facet_to_mask = fill(false,length(facet_to_cells))
   _fill_ghost_skeleton_mask!(facet_to_mask,facet_to_cells,cell_to_inoutcut,in_or_out.in_or_out)
 
   SkeletonTriangulation(model,facet_to_mask)
 end
 
-function _fill_ghost_skeleton_mask!(facet_to_mask,facet_to_cells::Table,cell_to_inoutcut,in_or_out)
-
-  nfacets = length(facet_to_cells)
-  for facet in 1:nfacets
+function _fill_ghost_skeleton_mask!(
+  facet_to_mask,facet_to_cells::Table,cell_to_inoutcut,in_or_out
+)
+  for facet in eachindex(facet_to_cells)
     a = facet_to_cells.ptrs[facet]
     b = facet_to_cells.ptrs[facet+1]
     ncells_around = b-a
@@ -539,5 +569,4 @@ function _fill_ghost_skeleton_mask!(facet_to_mask,facet_to_cells::Table,cell_to_
       facet_to_mask[facet] = true
     end
   end
-
 end

test/InterfacesTests/issue_115.jl

@@ -0,0 +1,54 @@
+# Tests for the fix to issue #115:
+#   https://github.com/gridap/GridapEmbedded.jl/issues/115
+#
+# When a compound CSG geometry (e.g. setdiff) is used, the conservative cell-level
+# classification can mark background cells as CUT even when no physical subcell of
+# the target domain lies in them. Before the fix, Triangulation(ACTIVE) and
+# GhostSkeleton would include these phantom cells, producing isolated elements.
+
+module Issue115Tests
+
+using Test
+using Gridap
+using Gridap.Geometry
+using Gridap.Adaptivity
+using GridapEmbedded
+using GridapEmbedded.LevelSetCutters
+
+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)
+
+# Two concentric circles: φ1 (radius r1=0.3), φ2 (radius r2 < r1 → empty setdiff)
+r1 = 0.3
+r2 = r1 - 0.5/n  # r2 < r1: Ω2 = setdiff(inside φ2, inside φ1) is geometrically empty
+φ1((x,y)) = (x-0.5)^2 + (y-0.5)^2 - r1^2
+φ2((x,y)) = (x-0.5)^2 + (y-0.5)^2 - r2^2
+
+V_φ = TestFESpace(model, ReferenceFE(lagrangian,Float64,1))
+φh1 = interpolate(φ1, V_φ)
+φh2 = interpolate(φ2, V_φ)
+
+geo1   = DiscreteGeometry(φh1, model, name="geo1")
+geo2   = DiscreteGeometry(φh2, model, name="geo2")
+Ω1_geo = setdiff(!geo1, geo2, name="Ω1")  # ring outside φ2 and inside φ1 (non-empty)
+Ω2_geo = setdiff(geo2,  geo1, name="Ω2")  # inside φ2 and outside φ1 (empty, r2 < r1)
+cutgeo = cut(model, union(Ω1_geo, Ω2_geo))
+
+# Empty domain: ACTIVE triangulation must contain no cells
+Ω2_act = Triangulation(cutgeo, ACTIVE, "Ω2")
+@test num_cells(Ω2_act) == 0
+
+# Empty domain: GhostSkeleton must contain no facets
+Λ2 = GhostSkeleton(cutgeo, ACTIVE_IN, "Ω2")
+@test num_cells(Λ2) == 0
+
+# Non-empty domain Ω1 must be unaffected
+Ω1_act = Triangulation(cutgeo, ACTIVE, "Ω1")
+@test num_cells(Ω1_act) > 0
+
+Λ1 = GhostSkeleton(cutgeo, ACTIVE_IN, "Ω1")
+@test num_cells(Λ1) > 0
+
+end # module

test/InterfacesTests/runtests.jl

@@ -12,4 +12,8 @@ using Test
 
 @testset "Cutters" begin include("CuttersTests.jl") end
 
+@testset "Issue Reproducers" begin
+  @testset "Issue115" begin include("issue_115.jl") end
+end
+
 end # module