Fix GridapEmbedded#115

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,10 +27,11 @@ AbstractTrees = "0.3.3, 0.4"
 Algoim = "0.2.2"
 Combinatorics = "1"
 CxxWrap = "0.16"
+DataStructures = "0.19"
 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"

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)

test/GridapEmbeddedTests/issues/issue_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),"Omega")
+writevtk(Triangulation(cutgeo,"Ω2"),"Omega2")
+writevtk(Triangulation(cutgeo,ACTIVE,"Ω2"),"Omega2_act")
+writevtk(Λ_Ω2,"Omega2_ghost_skel")

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