From 9e579875b4817ecbf8286677fe5da7188634a248 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Wed, 16 Oct 2024 12:51:46 +1100
Subject: [PATCH] Ghost face filtering

---
 scripts/Embedded/MWEs/jordi_embedded.jl |  6 +--
 scripts/Embedded/SubFacetSkeletons.jl   | 57 ++++++++++++++++++++++++-
 2 files changed, 58 insertions(+), 5 deletions(-)

diff --git a/scripts/Embedded/MWEs/jordi_embedded.jl b/scripts/Embedded/MWEs/jordi_embedded.jl
index 4071b651..1ecf53d6 100644
--- a/scripts/Embedded/MWEs/jordi_embedded.jl
+++ b/scripts/Embedded/MWEs/jordi_embedded.jl
@@ -9,7 +9,7 @@ import Gridap.Geometry: get_node_coordinates, collect1d
 include("../differentiable_trians.jl")
 
 order = 1
-n = 2
+n = 3
 N = 8
 
 # _model = CartesianDiscreteModel((0,1,0,1),(n,n))
@@ -33,7 +33,7 @@ V_φ = TestFESpace(model,reffe_scalar)
 # φh = interpolate(x->tan(-pi/4)*(x[1]-0.5)+(x[2]-0.5),V_φ) # Angled interface
 
 # φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.303,V_φ) # Circle
-φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2+(x[3]-0.5)^2)-0.303,V_φ) # Sphere
+φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2+(x[3]-0.5)^2)-0.25,V_φ) # Sphere
 
 fh = interpolate(x->cos(x[1]*x[2]),V_φ)
 
@@ -106,7 +106,7 @@ J2(φ) = ∫(fh)dΓ_AD
 dJ2_AD = gradient(J2,φh)
 dj2_AD = assemble_vector(dJ2_AD,V_φ) # TODO: Why is this +ve but AD on volume is -ve???
 
-dj2 - dj2_AD
+norm(dj2 - dj2_AD)
 
 ############################################################################################
 
diff --git a/scripts/Embedded/SubFacetSkeletons.jl b/scripts/Embedded/SubFacetSkeletons.jl
index 92225497..de9aa0a4 100644
--- a/scripts/Embedded/SubFacetSkeletons.jl
+++ b/scripts/Embedded/SubFacetSkeletons.jl
@@ -82,6 +82,57 @@ function generate_ghost_trian(
   return SkeletonTriangulation(plus,minus)
 end
 
+"""
+    get_ghost_mask(
+      face_trian::SubFacetTriangulation{Df,Dc},
+      face_model::UnstructuredDiscreteModel{Df,Dc}
+    ) where {Df,Dc}
+
+    get_ghost_mask(face_trian::SubFacetTriangulation) = get_ghost_mask(face_trian,get_active_model(face_trian))
+
+  Returns a mask for ghost faces. We define ghost faces as the interfaces between two 
+  different cut facets that are located in different background cells.
+
+  The second condition is important: In 3D, some cuts subcells may not be simplices. 
+  In this case, we simplexify the subcell. This creates extra cut interfaces that are 
+  interior to a background cell. These are not considered ghost faces.
+
+  - In 2D: Dc = 2, Df = 1 -> Ghost faces have dimension 0 (i.e interface points)
+  - In 3D: Dc = 3, Df = 2 -> Ghost faces have dimension 1 (i.e interface edges)
+"""
+function get_ghost_mask(
+  face_trian::SubFacetTriangulation{Df,Dc},
+  face_model::UnstructuredDiscreteModel{Df,Dc}
+) where {Df,Dc}
+  topo = get_grid_topology(face_model)
+  face_to_facets = get_faces(topo,Df-1,Df)
+
+  subfacets = face_trian.subfacets
+  facet_to_bgcell = subfacets.facet_to_bgcell
+  
+  n_faces = num_faces(topo,Df-1)
+  face_is_ghost = zeros(Bool,n_faces)
+  for face in 1:n_faces
+    facets = view(face_to_facets,face)
+    is_boundary = isone(length(facets))
+    if !is_boundary
+      @assert length(facets) == 2
+      bgcells = view(facet_to_bgcell,facets)
+      is_ghost = (bgcells[1] != bgcells[2])
+      face_is_ghost[face] = is_ghost
+    end
+  end
+
+  return face_is_ghost
+end
+
+function get_ghost_mask(
+  face_trian::SubFacetTriangulation
+)
+  face_model = get_active_model(face_trian)
+  return get_ghost_mask(face_trian,face_model)
+end
+
 struct SubFacetSkeletonTriangulation{Di,Df,Dp} <: Triangulation{Di,Dp}
   cell_skeleton::CompositeTriangulation{Di,Dp} # Interface -> BG Cell pair
   face_skeleton::SkeletonTriangulation{Di,Dp}  # Interface -> Cut Facet pair
@@ -90,10 +141,12 @@ struct SubFacetSkeletonTriangulation{Di,Df,Dp} <: Triangulation{Di,Dp}
   face_model::UnstructuredDiscreteModel{Df,Dp}
 end
 
-function SkeletonTriangulation(face_trian::SubFacetTriangulation)
+function Geometry.SkeletonTriangulation(face_trian::SubFacetTriangulation)
   bgmodel = get_background_model(face_trian)
   face_model = get_active_model(face_trian)
-  face_skeleton = SkeletonTriangulation(face_model)
+
+  ghost_mask = get_ghost_mask(face_trian,face_model)
+  face_skeleton = SkeletonTriangulation(face_model,ghost_mask)
   cell_skeleton = CompositeTriangulation(face_trian,face_skeleton)
   ghost_skeleton = generate_ghost_trian(cell_skeleton,bgmodel)
   return SubFacetSkeletonTriangulation(cell_skeleton,face_skeleton,ghost_skeleton,face_trian,face_model)