MWE fully working now

Project.toml

@@ -10,6 +10,7 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 CircularArrays = "7a955b69-7140-5f4e-a0ed-f168c5e2e749"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapDistributed = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
 GridapEmbedded = "8838a6a3-0006-4405-b874-385995508d5d"

scripts/Embedded/Examples/2d_thermal.jl

@@ -1,4 +1,4 @@
-using Gridap,GridapTopOpt
+using Gridap,GridapTopOpt, GridapSolvers
 using Gridap.Adaptivity, Gridap.Geometry
 using GridapEmbedded, GridapEmbedded.LevelSetCutters
 
@@ -10,9 +10,9 @@ mkpath(path)
 n = 51
 order = 1
 γ = 0.1
-max_steps = floor(Int,order*minimum(n)/10)
+max_steps = floor(Int,order*n/5)
 vf = 0.4
-α_coeff = 4max_steps*γ
+α_coeff = 3#4max_steps*γ
 iter_mod = 1
 
 _model = CartesianDiscreteModel((0,1,0,1),(n,n))
@@ -41,7 +41,7 @@ U_reg = TrialFESpace(V_reg,0)
 V_φ = TestFESpace(model,reffe_scalar)
 
 ## Levet-set function
-φh = interpolate(x->-cos(4π*x[1])*cos(4π*x[2])-0.2,V_φ)
+φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
   Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
@@ -71,11 +71,14 @@ Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ
 dVol(q,u,φ) = ∫(-1/vol_D*q/(norm ∘ (∇(φ))))Ωs.dΓ
 
 ## Setup solver and FE operators
+Pl = JacobiLinearSolver()
+ls = CGSolver(Pl;maxiter=1000,verbose=1,name="CG")
+
 state_collection = EmbeddedCollection(model,φh) do _
   # update_collection!(Ωs,φh)
   V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,0.0)
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
   (; 
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
@@ -85,28 +88,36 @@ end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method
-evo = CutFEMEvolve(V_φ,Ωs,dΩ,h)
+evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps)
 reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0))#InteriorPenalty(V_φ))
 ls_evo = UnfittedFEEvolution(evo,reinit)
 reinit!(ls_evo,φh)
 
 ## Hilbertian extension-regularisation problems
-## α = α_coeff*(h_refine/order)^2
-## a_hilb(p,q) =∫(α*∇(p)⋅∇(q) + p*q)dΩ;
-α = α_coeff*h_refine
-a_hilb(p,q) = ∫(α^2*∇(p)⋅∇(q) + p*q)dΩ;
+α = α_coeff*(h_refine/order)^2
+a_hilb(p,q) =∫(α*∇(p)⋅∇(q) + p*q)dΩ;
 vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
 
 ## Optimiser
+converged(m) = GridapTopOpt.default_al_converged(
+  m;
+  L_tol = 0.01*h_refine,
+  C_tol = 0.01
+)
+has_oscillations(m,os_it) = GridapTopOpt.default_has_oscillations(m,os_it;itlength=50,
+    itstart=100)
 optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;debug=true,
-  γ,verbose=true,constraint_names=[:Vol])
+  γ,verbose=true,constraint_names=[:Vol],converged)#,has_oscillations,reinit_mod=5)
 for (it,uh,φh,state) in optimiser
+  x_φ = get_free_dof_values(φh)
+  idx = findall(isapprox(0.0;atol=10^-10),x_φ)
+  !isempty(idx) && @warn "Boundary intersects nodes!"
   if iszero(it % iter_mod)
     writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"velh"=>FEFunction(V_φ,state.vel)])
-    writevtk(Ωs.Γ,path*"Gamma_out$it")
+    writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh])
   end
   write_history(path*"/history.txt",optimiser.history)
 end
 it = get_history(optimiser).niter; uh = get_state(pcfs)
 writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
-writevtk(Ωs.Γ,path*"Gamma_out$it")
+writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh])

src/GridapExtensions.jl

@@ -100,47 +100,49 @@ function Arrays.evaluate!(
   Polynomials._evaluate!(cache,fg,x,Val(false))
 end
 
-# function FESpaces._compute_cell_ids(uh,ttrian)
-#   strian = get_triangulation(uh)
-#   if strian === ttrian
-#     return collect(IdentityVector(Int32(num_cells(strian))))
-#   end
-#   @check is_change_possible(strian,ttrian)
-#   D = num_cell_dims(strian)
-#   sglue = get_glue(strian,Val(D))
-#   tglue = get_glue(ttrian,Val(D))
-#   @notimplementedif !isa(sglue,FaceToFaceGlue)
-#   @notimplementedif !isa(tglue,FaceToFaceGlue)
-#   scells = IdentityVector(Int32(num_cells(strian)))
-#   mcells = extend(scells,sglue.mface_to_tface)
-#   tcells = lazy_map(Reindex(mcells),tglue.tface_to_mface)
-#   # <-- Remove collect to keep PosNegReindex
-#   tcells = collect(tcells)
-#   return tcells
-# end
+# TODO: Below is dangerous, as it may break other Gridap methods,
+#   it is neccessary for now - see thermal_2d.jl problem
+function FESpaces._compute_cell_ids(uh,ttrian)
+  strian = get_triangulation(uh)
+  if strian === ttrian
+    return collect(IdentityVector(Int32(num_cells(strian))))
+  end
+  @check is_change_possible(strian,ttrian)
+  D = num_cell_dims(strian)
+  sglue = get_glue(strian,Val(D))
+  tglue = get_glue(ttrian,Val(D))
+  @notimplementedif !isa(sglue,FaceToFaceGlue)
+  @notimplementedif !isa(tglue,FaceToFaceGlue)
+  scells = IdentityVector(Int32(num_cells(strian)))
+  mcells = extend(scells,sglue.mface_to_tface)
+  tcells = lazy_map(Reindex(mcells),tglue.tface_to_mface)
+  # <-- Remove collect to keep PosNegReindex
+  # tcells = collect(tcells)
+  return tcells
+end
 
 # New dispatching
-# function Arrays.lazy_map(k::Reindex,ids::Arrays.LazyArray{<:Fill{<:PosNegReindex}})
-#   k_posneg = ids.maps.value
-#   posneg_partition = ids.args[1]
-#   pos_values = lazy_map(Reindex(k.values),k_posneg.values_pos)
-#   pos_values, neg_values = Geometry.pos_neg_data(pos_values,posneg_partition)
-#   # println("Byee ::: $(eltype(pos_values)) --- $(eltype(neg_values))")
-#   lazy_map(PosNegReindex(pos_values,neg_values),posneg_partition)
-# end
-
-# function Arrays.lazy_map(k::Reindex,ids::Arrays.AppendedArray)
-#   a = lazy_map(k,ids.a)
-#   b = lazy_map(k,ids.b)
-#   # println("Hello ::: $(eltype(a)) --- $(eltype(b))")
-#   return lazy_append(a,b)
-# end
-
-# using ForwardDiff
-
-# function Arrays.evaluate!(result,k::AutoDiffMap,ydual,x,cfg::ForwardDiff.GradientConfig{T}) where T
-#   @notimplementedif ForwardDiff.chunksize(cfg) != length(x)
-#   @notimplementedif length(result) != length(x)
-#   !isempty(x) && ForwardDiff.extract_gradient!(T, result, ydual) # <-- Watch for empty cell contributions
-#   return result
-# end
+function Arrays.lazy_map(k::Reindex,ids::Arrays.LazyArray{<:Fill{<:PosNegReindex}})
+  k_posneg = ids.maps.value
+  posneg_partition = ids.args[1]
+  pos_values = lazy_map(Reindex(k.values),k_posneg.values_pos)
+  pos_values, neg_values = Geometry.pos_neg_data(pos_values,posneg_partition)
+  # println("Byee ::: $(eltype(pos_values)) --- $(eltype(neg_values))")
+  lazy_map(PosNegReindex(pos_values,neg_values),posneg_partition)
+end
+
+function Arrays.lazy_map(k::Reindex,ids::Arrays.AppendedArray)
+  a = lazy_map(k,ids.a)
+  b = lazy_map(k,ids.b)
+  # println("Hello ::: $(eltype(a)) --- $(eltype(b))")
+  return lazy_append(a,b)
+end
+
+using ForwardDiff
+
+function Arrays.evaluate!(result,k::AutoDiffMap,ydual,x,cfg::ForwardDiff.GradientConfig{T}) where T
+  @notimplementedif ForwardDiff.chunksize(cfg) != length(x)
+  @notimplementedif length(result) != length(x)
+  !isempty(x) && ForwardDiff.extract_gradient!(T, result, ydual) # <-- Watch for empty cell contributions
+  return result
+end