diff --git a/Project.toml b/Project.toml
index f277a42..63c241b 100755
--- a/Project.toml
+++ b/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"
diff --git a/scripts/Embedded/Examples/2d_thermal copy.jl b/scripts/Embedded/Examples/2d_thermal copy.jl
deleted file mode 100644
index 5938f84..0000000
--- a/scripts/Embedded/Examples/2d_thermal copy.jl	
+++ /dev/null
@@ -1,142 +0,0 @@
-using Gridap,GridapTopOpt
-using Gridap.Adaptivity, Gridap.Geometry
-using GridapEmbedded, GridapEmbedded.LevelSetCutters
-
-using GridapTopOpt: StateParamIntegrandWithMeasure
-
-path="./results/UnfittedFEM_thermal_compliance_ALM/"
-n = 51
-order = 1
-γ = 0.1
-γ_reinit = 0.5
-max_steps = floor(Int,order*minimum(n)/10)
-tol = 1/(5*order^2)/minimum(n)
-vf = 0.4
-α_coeff = 4max_steps*γ
-iter_mod = 1
-
-_model = CartesianDiscreteModel((0,1,0,1),(n,n))
-base_model = UnstructuredDiscreteModel(_model)
-ref_model = refine(base_model, refinement_method = "barycentric")
-model = ref_model.model
-el_Δ = get_el_Δ(_model)
-h = maximum(el_Δ)
-f_Γ_D(x) = (x[1] ≈ 0.0 && (x[2] <= 0.2 + eps() || x[2] >= 0.8 - eps()))
-f_Γ_N(x) = (x[1] ≈ 1 && 0.4 - eps() <= x[2] <= 0.6 + eps())
-update_labels!(1,model,f_Γ_D,"Gamma_D")
-update_labels!(2,model,f_Γ_N,"Gamma_N")
-
-## Triangulations and measures
-Ω = Triangulation(model)
-Γ_N = BoundaryTriangulation(model,tags="Gamma_N")
-dΩ = Measure(Ω,2*order)
-dΓ_N = Measure(Γ_N,2*order)
-vol_D = sum(∫(1)dΩ)
-
-## Levet-set function space and derivative regularisation space
-reffe_scalar = ReferenceFE(lagrangian,Float64,order)
-V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"])
-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_φ)
-geo = DiscreteGeometry(φh,model)
-cutgeo = cut(model,geo)
-Ωs = EmbeddedCollection(model,φh) do cutgeo
-  Ωin = Triangulation(cutgeo,PHYSICAL)
-  Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
-  (; 
-    :Ωin  => Ωin,
-    :dΩin => Measure(Ωin,2*order),
-    :Γg   => Γg,
-    :dΓg  => Measure(Γg,2*order),
-    :n_Γg => get_normal_vector(Γg)  
-  )
-end  
-
-## Weak form
-const γg = 0.1
-a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + ∫((γg*h)*jump(Ωs.n_Γg⋅∇(v))*jump(Ωs.n_Γg⋅∇(u)))Ωs.dΓg
-l(v,φ) = ∫(v)dΓ_N
-
-## Optimisation functionals
-J(u,φ) = a(u,u,φ)
-Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ
-
-## Setup solver and FE operators
-state_collection = EmbeddedCollection(model,φh) do cutgeo
-  Ωact = Triangulation(cutgeo,ACTIVE)
-  V = TestFESpace(Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
-  U = TrialFESpace(V,0.0)
-
-  Ωin = Triangulation(cutgeo,PHYSICAL)
-  Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
-  dΩin = Measure(Ωin,2*order)
-  dΓg  = Measure(Γg,2*order)
-  n_Γg = get_normal_vector(Γg)  
-
-  dΩact = Measure(Ωact,2)
-  Ωact_out = Triangulation(cutgeo,ACTIVE_OUT)
-  dΩact_out = Measure(Ωact_out,2)
-
-  a(u,v,φ) = ∫(∇(v)⋅∇(u))dΩin + ∫((γg*h)*jump(n_Γg⋅∇(v))*jump(n_Γg⋅∇(u)))dΓg
-  l( v,φ) = ∫(v)dΓ_N
-  J(u,φ) = a(u,u,φ)
-  Vol(u,φ) = ∫(1/vol_D)dΩin - ∫(vf/vol_D)dΩact - ∫(vf/vol_D)dΩact_out
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
-  (; 
-    :state_map => state_map,
-    :J => StateParamIntegrandWithMeasure(J,state_map),
-    :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),(Vol,))
-  )
-end  
-pcfs = EmbeddedPDEConstrainedFunctionals(state_collection)
-
-evaluate!(pcfs,φh)
-
-using Gridap.Arrays
-uh = u
-ttrian = Ω
-strian = get_triangulation(uh)
-
-D = num_cell_dims(strian)
-sglue = get_glue(strian,Val(D))
-tglue = get_glue(ttrian,Val(D))
-
-scells = Arrays.IdentityVector(Int32(num_cells(strian)))
-mcells = extend(scells,sglue.mface_to_tface)
-tcells = lazy_map(Reindex(mcells),tglue.tface_to_mface)
-collect(tcells)
-
-
-## Evolution Method
-evo = CutFEMEvolve(V_φ,Ωs,dΩ,h)
-reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=InteriorPenalty(V_φ))
-ls_evo = UnfittedFEEvolution(evo,reinit)
-reinit!(ls_evo,φh)
-
-## Hilbertian extension-regularisation problems
-α = α_coeff*maximum(el_Δ)
-a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ;
-vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
-
-Ωs(φh)
-
-evaluate!(pcfs,φh)
- 
-# ## Optimiser
-# optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
-#   γ,γ_reinit,verbose=true,constraint_names=[:Vol])
-# for (it,uh,φh) in optimiser
-#   if iszero(it % iter_mod)
-#     writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
-#     writevtk(Ωs.Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)])
-#   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",cellfields=["normal"=>get_normal_vector(_Γ)])
\ No newline at end of file
diff --git a/scripts/Embedded/Examples/2d_thermal.jl b/scripts/Embedded/Examples/2d_thermal.jl
index 8c86ae1..cbb112f 100644
--- a/scripts/Embedded/Examples/2d_thermal.jl
+++ b/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")
\ No newline at end of file
+writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh])
\ No newline at end of file
diff --git a/src/GridapExtensions.jl b/src/GridapExtensions.jl
index bf31e29..db7abc2 100644
--- a/src/GridapExtensions.jl
+++ b/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
\ No newline at end of file
+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
\ No newline at end of file