Refactored stuff

scripts/Embedded/MWEs/jordi_embedded.jl

@@ -6,13 +6,11 @@ using GridapEmbedded.LevelSetCutters
 using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Adaptivity
 import Gridap.Geometry: get_node_coordinates, collect1d
 
-include("../differentiable_trians.jl")
-
 order = 1
 n = 15
 N = 8
 
-# _model = CartesianDiscreteModel((0,1,0,1),(n,n))
+#_model = CartesianDiscreteModel((0,1,0,1),(n,n))
 _model = CartesianDiscreteModel((0,1,0,1,0,1),(n,n,n))
 cd = Gridap.Geometry.get_cartesian_descriptor(_model)
 base_model = UnstructuredDiscreteModel(_model)
@@ -31,16 +29,16 @@ V_φ = TestFESpace(model,reffe_scalar)
 # φh = interpolate(x->(1+0.25)abs(x[1]-0.5)+0abs(x[2]-0.5)-0.25,V_φ) # Straight lines with scaling
 # φh = interpolate(x->abs(x[1]-0.5)+0abs(x[2]-0.5)-0.25/(1+0.25),V_φ) # Straight lines without scaling
 # φ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.5223,V_φ) # Circle
+#φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.5223,V_φ) # Circle
 
 # φh = interpolate(x->max(abs(x[1]-0.5),abs(x[2]-0.5),abs(x[3]-0.5))-0.25,V_φ) # Square prism
 # φh = interpolate(x->((x[1]-0.5)^N+(x[2]-0.5)^N+(x[3]-0.5)^N)^(1/N)-0.25,V_φ) # Curved corner example
 # φh = interpolate(x->abs(x[1]-0.5)+abs(x[2]-0.5)+abs(x[3]-0.5)-0.25-0/n/10,V_φ) # Diamond prism
 # φh = interpolate(x->(1+0.25)abs(x[1]-0.5)+0abs(x[2]-0.5)+0abs(x[3]-0.5)-0.25,V_φ) # Straight lines with scaling
 # φh = interpolate(x->abs(x[1]-0.5)+0abs(x[2]-0.5)+0abs(x[3]-0.5)-0.25/(1+0.25),V_φ) # Straight lines without scaling
 # φh = interpolate(x->tan(-pi/3)*(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+(x[3]-0.5)^2)-0.53,V_φ) # Sphere
-φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])*cos(2π*x[3])-0.11,V_φ) # "Regular" LSF
+φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2+(x[3]-0.5)^2)-0.53,V_φ) # Sphere
+#φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])*cos(2π*x[3])-0.11,V_φ) # "Regular" LSF
 
 fh = interpolate(x->cos(x[1]*x[2]),V_φ)
 
@@ -84,7 +82,8 @@ println("dj1 error = $(norm(dj_vec_out-dj_exp_vec))")
 
 ############################################################################################
 
-include("../SubFacetSkeletons.jl")
+using GridapTopOpt: get_subfacet_normal_vector, get_ghost_normal_vector
+using GridapTopOpt: get_conormal_vector, get_tangent_vector
 
 Γ = EmbeddedBoundary(cutgeo)
 Λ = Skeleton(Γ)
@@ -103,8 +102,8 @@ m_k = get_conormal_vector(Λ)
 
 fh = interpolate(x->1.0,V_φ)
 ∇ˢφ = Operation(abs)(n_S ⋅ ∇(φh).plus)
-_n = get_normal_vector(Γ)
-dJ2(w) = ∫((-_n⋅∇(fh))*w/(norm ∘ (∇(φh))))dΓ + ∫(-n_S ⋅ (jump(fh*m_k) * mean(w) / ∇ˢφ))dΛ
+n_Γ = get_normal_vector(Γ)
+dJ2(w) = ∫((-n_Γ⋅∇(fh))*w/(norm ∘ (∇(φh))))dΓ + ∫(-n_S ⋅ (jump(fh*m_k) * mean(w) / ∇ˢφ))dΛ
 dj2 = assemble_vector(dJ2,V_φ)
 
 diff_Γ = DifferentiableTriangulation(Γ)
@@ -125,7 +124,7 @@ m_k_Σ = get_conormal_vector(Σ)
 ∇ˢφ_Σ = Operation(abs)(n_S_Σ ⋅ ∇(φh))
 
 dΣ = Measure(Σ,2*order)
-dJ3(w) = dJ2(w) + ∫(-n_S_Σ ⋅ (fh*m_k_Σ * w / ∇ˢφ_Σ))dΣ
+dJ3(w) = ∫(-n_S_Σ ⋅ (fh*m_k_Σ * w / ∇ˢφ_Σ))dΣ + dJ2(w)
 dj3 = assemble_vector(dJ3,V_φ)
 
 println("dj2 error = $(norm(dj3 - dj2_AD))")
@@ -171,4 +170,26 @@ writevtk(
 writevtk(
   Γ,"results/gamma";
   append=false
-)
+)
+
+Boundary(Γ)
+
+face_trian = Γ
+
+bgmodel = get_background_model(face_trian)
+face_model = get_active_model(face_trian)
+
+face_boundary  = BoundaryTriangulation(face_model)
+cell_boundary  = Geometry.CompositeTriangulation(face_trian,face_boundary)
+ghost_boundary = GridapTopOpt.generate_ghost_trian(cell_boundary,bgmodel)
+interface_sign = get_interface_sign(face_trian,cell_boundary,ghost_boundary)
+
+GridapTopOpt.get_ghost_facet_normal(cell_boundary,ghost_boundary)
+GridapTopOpt.get_subfacet_facet_normal(cell_boundary,face_trian)
+s = GridapTopOpt.get_interface_sign(cell_boundary,face_trian,ghost_boundary)
+
+Λbg = Skeleton(bgmodel)
+
+nbg = get_normal_vector(Λbg)
+
+

scripts/Embedded/SubFacetSkeletons.jl

@@ -237,18 +237,12 @@ function CellData.change_domain(
 end
 
 # Normal vector to the cut facets , n_∂Ω
-function get_subfacet_normal_vector(trian::SubFacetSkeletonTriangulation{0})
+function get_subfacet_normal_vector(trian::SubFacetSkeletonTriangulation)
   n_∂Ω = get_normal_vector(trian.face_trian)
   plus = change_domain(n_∂Ω.plus,trian,ReferenceDomain())
   minus = change_domain(n_∂Ω.minus,trian,ReferenceDomain())
   return SkeletonPair(plus,minus)
 end
-function get_subfacet_normal_vector(trian::SubFacetSkeletonTriangulation{1})
-  n_∂Ω = get_normal_vector(trian.face_trian)
-  plus = change_domain(n_∂Ω.plus,trian,ReferenceDomain())
-  minus = change_domain(n_∂Ω.minus,trian,ReferenceDomain()) # This needs to be postive for mk to be correctly oriented in 3d
-  return SkeletonPair(plus,minus)
-end
 
 # Normal vector to the ghost facets, n_k
 # function get_ghost_normal_vector(trian::SubFacetSkeletonTriangulation)