2d thermal MWE

scripts/Embedded/Bugs/autodiff_u.jl

@@ -0,0 +1,30 @@
+using Gridap
+using Gridap.CellData
+
+function lazy_collect(a)
+    c = array_cache(a)
+    for i in eachindex(a)
+        getindex!(c,a,i)
+    end
+end
+
+model = CartesianDiscreteModel((0,1,0,1),(2,2))
+
+Ω_space = Triangulation(model,[1,2])
+#Ω_space = Triangulation(model,[3,4])
+
+reffe = ReferenceFE(lagrangian,Float64,1)
+V = FESpace(Ω_space,reffe)
+assem = SparseMatrixAssembler(V,V)
+
+Ω = Triangulation(model)
+dΩ = Measure(Ω,2)
+j(u) = ∫(1.0)dΩ
+
+uh = zero(V)
+c = get_contribution(gradient(j,uh),Ω)
+
+collect(c) # Fails
+lazy_collect(c)
+
+assemble_vector(assem,collect_cell_vector(V,gradient(j,uh)))

scripts/Embedded/Examples/2d_thermal.jl

@@ -6,13 +6,13 @@ using GridapTopOpt: StateParamIntegrandWithMeasure
 
 path="./results/UnfittedFEM_thermal_compliance_ALM/"
 mkpath(path)
-n = 51
+n = 101
 order = 1
 γ = 0.1
 γ_reinit = 0.5
 max_steps = floor(Int,order*minimum(n)/10)
 vf = 0.4
-α_coeff = 2max_steps*γ
+α_coeff = 4max_steps*γ
 iter_mod = 1
 
 _model = CartesianDiscreteModel((0,1,0,1),(n,n))
@@ -44,17 +44,19 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(4π*x[1])*cos(4π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γg = DifferentiableTriangulation(GhostSkeleton(cutgeo))
-  n_Γg = get_normal_vector(Γg)
   Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γg = GhostSkeleton(cutgeo)
+  n_Γg = get_normal_vector(Γg)
+  Ωact = Triangulation(cutgeo,ACTIVE)
   (; 
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
     :dΓg  => Measure(Γg,2*order),
     :n_Γg => get_normal_vector(Γg),
     :Γ    => Γ,
-    :dΓ   => Measure(Γ,2*order)
+    :dΓ   => Measure(Γ,2*order),
+    :Ωact => Ωact
   )
 end  
 
@@ -64,14 +66,14 @@ a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + ∫((γg*h)*jump(Ωs.n_Γg⋅∇(v))
 l(v,φ) = ∫(v)dΓ_N
 
 ## Optimisation functionals
-J(u,φ) = ∫(∇(u)⋅∇(u))Ωs.dΩin
+J(u,φ) = a(u,u,φ)
 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
-state_collection = EmbeddedCollection(model,φh) do cutgeo
-  Ωact = Triangulation(cutgeo,ACTIVE)
-  V = TestFESpace(Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
+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)
   (; 
@@ -84,32 +86,27 @@ pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method
 evo = CutFEMEvolve(V_φ,Ωs,dΩ,h)
-reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=InteriorPenalty(V_φ))
+reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0))#InteriorPenalty(V_φ))
 ls_evo = UnfittedFEEvolution(evo,reinit)
-# reinit!(ls_evo,φh)
+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/order)^2
+## a_hilb(p,q) =∫(α*∇(p)⋅∇(q) + p*q)dΩ;
 α = α_coeff*h_refine
 a_hilb(p,q) = ∫(α^2*∇(p)⋅∇(q) + p*q)dΩ;
 vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
 
-
-## Forward problem and derivatives
-_j,_c,_dJ,_dC = evaluate!(pcfs,φh)
-writevtk(Ω,path*"check_dJ",cellfields=["dJ"=>FEFunction(dJ,U_reg)])
-
 ## Optimiser
-optimiser = HilbertianProjection(pcfs,ls_evo,vel_ext,φh;debug=true,
+optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;debug=true,
   γ,γ_reinit,verbose=true,constraint_names=[:Vol])
 for (it,uh,φh,state) in optimiser
   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",cellfields=["normal"=>get_normal_vector(Ωs.Γ)])
+    writevtk(Ωs.Γ,path*"Gamma_out$it")
   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(Ωs.Γ)])
+writevtk(Ωs.Γ,path*"Gamma_out$it")

src/Embedded/DifferentiableTriangulations.jl

@@ -45,17 +45,24 @@ end
 function FESpaces._change_argument(
   op,f,trian::DifferentiableTriangulation,uh::SingleFieldFEFunction
 )
+  update = matching_level_set(trian,uh)
   U = get_fe_space(uh)
   function g(cell_u)
     cf = CellField(U,cell_u)
-    update_trian!(trian,cf)
+    update ? update_trian!(trian,cf) : update_trian!(trian,nothing)
     cell_grad = f(cf)
     update_trian!(trian,nothing) # TODO: experimental
     get_contribution(cell_grad,trian)
   end
   g
 end
 
+function matching_level_set(trian::DifferentiableTriangulation,uh)
+  uh_trian = get_triangulation(uh)
+  bgmodel  = get_background_model(uh_trian)
+  return num_cells(uh_trian) == num_cells(bgmodel)
+end
+
 function FESpaces._compute_cell_ids(uh,ttrian::DifferentiableTriangulation)
   FESpaces._compute_cell_ids(uh,ttrian.trian)
 end
@@ -389,10 +396,11 @@ end
 function FESpaces._change_argument(
   op,f,trian::AppendedTriangulation,uh::SingleFieldFEFunction
 )
+  update = matching_level_set(trian,uh)
   U = get_fe_space(uh)
   function g(cell_u)
     cf = CellField(U,cell_u)
-    update_trian!(trian,cf)
+    update ? update_trian!(trian,cf) : update_trian!(trian,nothing)
     cell_grad = f(cf)
     update_trian!(trian,nothing) # TODO: experimental
     get_contribution(cell_grad,trian)
@@ -405,3 +413,13 @@ function FESpaces._compute_cell_ids(uh,ttrian::AppendedTriangulation)
   ids_b = FESpaces._compute_cell_ids(uh,ttrian.b)
   lazy_append(ids_a,ids_b)
 end
+
+function matching_level_set(trian::AppendedTriangulation,uh)
+  a = matching_level_set(trian.a,uh)
+  b = matching_level_set(trian.b,uh)
+  return a || b
+end
+
+function matching_level_set(trian::Triangulation,uh)
+  false
+end

src/GridapExtensions.jl

@@ -83,3 +83,64 @@ end
 function get_local_assembly_strategy(a::MultiField.BlockSparseMatrixAssembler)
   return get_local_assembly_strategy(first(a.block_assemblers))
 end
+
+# Fix for isbitstype bug in Gridap.Polynomials
+function Arrays.return_cache(
+  fg::Fields.FieldGradientArray{1,Polynomials.MonomialBasis{D,V}},
+  x::AbstractVector{<:Point}) where {D,V}
+  xi = testitem(x)
+  T = gradient_type(V,xi)
+  Polynomials._return_cache(fg,x,T,Val(false))
+end
+
+function Arrays.evaluate!(
+  cache,
+  fg::Fields.FieldGradientArray{1,Polynomials.MonomialBasis{D,V}},
+  x::AbstractVector{<:Point}) where {D,V}
+  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
+
+# 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

src/GridapTopOpt.jl

@@ -12,7 +12,7 @@ using DelimitedFiles, Printf
 using Gridap
 using Gridap.Helpers, Gridap.Algebra, Gridap.TensorValues
 using Gridap.Geometry, Gridap.CellData, Gridap.Fields, Gridap.Arrays
-using Gridap.ReferenceFEs, Gridap.FESpaces,  Gridap.MultiField
+using Gridap.ReferenceFEs, Gridap.FESpaces,  Gridap.MultiField, Gridap.Polynomials
 
 using Gridap.Geometry: get_faces, num_nodes
 using Gridap.FESpaces: get_assembly_strategy

src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl

@@ -10,7 +10,7 @@ Burman et al. (2017). DOI: `10.1016/j.cma.2017.09.005`.
 - `dΩ_bg::C`: Measure for integration
 - `space::B`: Level-set FE space
 - `assembler::Assembler`: FE assembler
-- `params::D`: Tuple of stabilisation parameter `Γg`, mesh size `h`, and
+- `params::D`: Tuple of stabilisation parameter `γg`, mesh size `h`, and
   max steps `max_steps`
 - `cache`: Cache for evolver, initially `nothing`.
 
@@ -30,29 +30,29 @@ mutable struct CutFEMEvolve{A,B,C} <: Evolver
   cache
   function CutFEMEvolve(V_φ::B,Ωs::EmbeddedCollection,dΩ_bg::A,h::Real;
       max_steps=10,
-      Γg = 0.1,
+      γg = 0.1,
       ode_ls = LUSolver(),
       ode_nl = NLSolver(ode_ls, show_trace=false, method=:newton, iterations=10),
       ode_solver = MutableRungeKutta(ode_nl, ode_ls, 0.1, :DIRK_CrankNicolson_2_2),
       assembler=SparseMatrixAssembler(V_φ,V_φ)) where {A,B}
-    if !(:F_act ∈ keys(Ωs.objects))
-      @warn "Expected triangulation ':F_act' not found in the 
-      EmbeddedCollection. This and the corresponding measure ':dF_act' have been
+    if !((:dΓg,:n_Γg) ⊆ keys(Ωs.objects))
+      @warn "Expected triangulation ':dΓg, :n_Γg' not found in the 
+      EmbeddedCollection. These and the corresponding triangulation ':Γg' have been
       added to the recipe list. 
       
-      Ensure that you are not using ':F_act' under a different
+      Ensure that you are not using ':Γg' under a different
       name to avoid additional computation for cutting."
       function F_act_recipe(cutgeo)
-        Ω_act = Triangulation(cutgeo,ACTIVE)
-        F_act = SkeletonTriangulation(Ω_act)
+        Γg = GhostSkeleton(cutgeo)
         (; 
-          :F_act => F_act,
-          :dF_act => Measure(F_act,2get_order(V_φ))
+          :Γg => Γg,
+          :dΓg => Measure(Γg,2get_order(V_φ)),
+          :n_Γg => get_normal_vector(Γg),
         )
       end
       add_recipe!(Ωs,F_act_recipe)
     end
-    params = (;Γg,h,max_steps)
+    params = (;γg,h,max_steps)
     new{A,B,typeof(params)}(ode_solver,Ωs,dΩ_bg,V_φ,assembler,params,nothing)
   end
 end
@@ -65,18 +65,28 @@ get_measure(s::CutFEMEvolve) = s.dΩ_bg
 get_params(s::CutFEMEvolve) = s.params
 get_cache(s::CutFEMEvolve) = s.cache
 
-function get_transient_operator(φh,velh,s::CutFEMEvolve)
+function get_transient_operator(φh,velh,s::CutFEMEvolve;use_full_skeleton=false)
   Ωs, V_φ, dΩ_bg, assembler, params = s.Ωs, s.space, s.dΩ_bg, s.assembler, s.params
-  Γg, h = params.Γg, params.h
+  γg, h = params.γg, params.h
   ϵ = 1e-20
 
-  update_collection!(Ωs,φh)
-  F_act = Ωs.F_act
-  dF_act = Ωs.dF_act
-  n = get_normal_vector(F_act)
+  # NOTE/TODO: The sparsity pattern of the Jacobian changes with the level-set function
+  #   because of the `∫(γg*h^2*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg` term, where
+  #   dΓg is the measure on the ghost skeleton. As such, the first time we compute the
+  #   operator we use the full background mesh skeleton so that the sparsity pattern
+  #   of the Jacobian is the worst possible. Subsequent iterations will re-use the sparsity
+  #   pattern but less integration is done. The trade-off here is memory cost vs. integration.
+  if use_full_skeleton
+    Γg = SkeletonTriangulation(get_triangulation(dΩ_bg.quad))
+    dΓg = Measure(Γg,2get_order(V_φ))
+    n_Γg = get_normal_vector(Γg)
+  else
+    update_collection!(Ωs,φh)
+    dΓg, n_Γg = Ωs.dΓg, Ωs.n_Γg
+  end
 
   β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh))
-  stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ_bg + ∫(Γg*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dF_act
+  stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ_bg + ∫(γg*h^2*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg
   mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ_bg
   forcing(t,v) = ∫(0v)dΩ_bg
   Ut_φ = TransientTrialFESpace(V_φ)
@@ -102,7 +112,7 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing)
   h, max_steps = params.h, params.max_steps
 
   # Setup FE operator and solver
-  ode_op = get_transient_operator(φh,velh,s)
+  ode_op = get_transient_operator(φh,velh,s;use_full_skeleton=true)
   dt = _compute_Δt(h,γ,get_free_dof_values(velh))
   ode_solver.dt = dt
   ode_sol = solve(ode_solver,ode_op,0.0,dt*max_steps,φh)

src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl

@@ -44,7 +44,7 @@ mutable struct StabilisedReinit{A,B,C,D} <: Reinitialiser
       Ensure that you are not using ':dΓ' under a different
       name to avoid additional computation for cutting."
       function dΓ_recipe(cutgeo)
-        Γ = EmbeddedBoundary(cutgeo)
+        Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
         (; 
           :dΓ => Measure(Γ,2get_order(V_φ))
         )
@@ -111,6 +111,7 @@ function get_residual_and_jacobian(s::StabilisedReinit{ArtificialViscosity},φh)
   b(w,v) = ∫((sign ∘ w)*v)dΩ_bg
   res(u,v) = a(u,u,v) - b(u,v)
   jac(u,du,v) = a(u,du,v)
+  # jac(u,du,v) = jacobian(res,[u,v],1)
   return res,jac
 end
 
@@ -137,6 +138,6 @@ function get_residual_and_jacobian(s::StabilisedReinit{InteriorPenalty},φh)
   b(w,v) = ∫((sign ∘ w)*v)dΩ_bg
   res(u,v) = a(u,u,v) - b(u,v)
   jac(u,du,v) = a(u,du,v)
-
+  # jac(u,du,v) = jacobian(res,[u,v],1)
   return res,jac
 end

test/seq/UnfittedEvolutionTests/CutFEMEvolveTest.jl

@@ -23,9 +23,13 @@ V_φ = TestFESpace(model,reffe_scalar)
 
 Ωs = EmbeddedCollection(model,φh) do cutgeo
   Γ = EmbeddedBoundary(cutgeo)
+  Γg = GhostSkeleton(cutgeo)
   (; 
     :Γ  => Γ,
     :dΓ => Measure(Γ,2*order),
+    :Γg => Γg,
+    :dΓg => Measure(Γg,2*order),
+    :n_Γg => get_normal_vector(Γg)
   )
 end