fcm approach

scripts/Embedded/Examples/FCM_2d_thermal.jl

@@ -0,0 +1,113 @@
+using Gridap,GridapTopOpt, GridapSolvers
+using Gridap.Adaptivity, Gridap.Geometry
+using GridapEmbedded, GridapEmbedded.LevelSetCutters
+
+using GridapTopOpt: StateParamIntegrandWithMeasure
+
+path="./results/FCM_thermal_compliance_ALM/"
+rm(path,force=true,recursive=true)
+mkpath(path)
+n = 50
+order = 1
+γ = 0.1
+max_steps = floor(Int,order*n/5)
+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_Δ)
+h_refine = maximum(el_Δ)/2
+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.4,V_φ)
+Ωs = EmbeddedCollection(model,φh) do cutgeo,_
+  Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_IN),V_φ)
+  Ωout = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ)
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
+  Γg = GhostSkeleton(cutgeo)
+  n_Γg = get_normal_vector(Γg)
+  Ωact = Triangulation(cutgeo,ACTIVE)
+  (; 
+    :Ωin  => Ωin,
+    :dΩin => Measure(Ωin,2*order),
+    :Ωout  => Ωout,
+    :dΩout => Measure(Ωout,2*order),
+    :Γg   => Γg,
+    :dΓg  => Measure(Γg,2*order),
+    :n_Γg => get_normal_vector(Γg),
+    :Γ    => Γ,
+    :dΓ   => Measure(Γ,2*order),
+    :Ωact => Ωact
+  )
+end
+
+## Weak form
+const ϵ = 1e-3
+a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + ∫(ϵ*∇(v)⋅∇(u))Ωs.dΩout
+l(v,φ) = ∫(v)dΓ_N
+
+## Optimisation functionals
+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
+V = TestFESpace(Ω,reffe_scalar;dirichlet_tags=["Gamma_D"])
+U = TrialFESpace(V,0.0)
+state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
+pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dC=(dVol,))
+
+## Evolution Method
+evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps)
+reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0))
+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Ω;
+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
+)
+optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;debug=true,
+  γ,verbose=true,constraint_names=[:Vol],converged)
+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.Ω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.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh])

src/Embedded/DifferentiableTriangulations.jl

@@ -11,58 +11,57 @@ methods to compute derivatives wrt deformations of the embedded mesh.
 To do so, it propagates dual numbers into the geometric maps mapping cut subcells/subfacets 
 to the background mesh.
 """
-mutable struct DifferentiableTriangulation{Dc,Dp,A} <: Triangulation{Dc,Dp}
+mutable struct DifferentiableTriangulation{Dc,Dp,A,B} <: Triangulation{Dc,Dp}
   trian :: A
+  fe_space :: B
   cell_values
   caches
   function DifferentiableTriangulation(
-    trian :: Triangulation{Dc,Dp},cell_values,caches
+    trian :: Triangulation{Dc,Dp},
+    fe_space :: FESpace,
+    cell_values,caches
   ) where {Dc,Dp}
     A = typeof(trian)
-    new{Dc,Dp,A}(trian,cell_values,caches)
+    B = typeof(fe_space)
+    new{Dc,Dp,A,B}(trian,fe_space,cell_values,caches)
   end
 end
 
 function DifferentiableTriangulation(
-  trian::Union{<:SubCellTriangulation,<:SubFacetTriangulation}
+  trian::Union{<:SubCellTriangulation,<:SubFacetTriangulation},
+  fe_space
 )
   caches = precompute_autodiff_caches(trian)
-  return DifferentiableTriangulation(trian,nothing,caches)
+  return DifferentiableTriangulation(trian,fe_space,nothing,caches)
 end
 
-(t::DifferentiableTriangulation)(φh) = update_trian!(t,φh)
+(t::DifferentiableTriangulation)(φh) = update_trian!(t,get_fe_space(φh),φh)
 
-function update_trian!(trian::DifferentiableTriangulation,φh)
+function update_trian!(trian::DifferentiableTriangulation,U::FESpace,φh)
+  ~(U === trian.fe_space) && return trian
   trian.cell_values = extract_dualized_cell_values(trian.trian,φh)
   return trian
 end
 
-function update_trian!(trian::DifferentiableTriangulation,::Nothing)
+function update_trian!(trian::DifferentiableTriangulation,::FESpace,::Nothing)
   trian.cell_values = nothing
   return trian
 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 ? update_trian!(trian,cf) : update_trian!(trian,nothing)
+    update_trian!(trian,U,cf)
     cell_grad = f(cf)
-    update_trian!(trian,nothing) # TODO: experimental
+    update_trian!(trian,U,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
@@ -375,34 +374,34 @@ end
 # We only need to propagate the dual numbers to the CUT cells, which is what the 
 # following implementation does:
 
-DifferentiableTriangulation(trian::Triangulation) = trian
+DifferentiableTriangulation(trian::Triangulation,fe_space) = trian
 
 function DifferentiableTriangulation(
-  trian::AppendedTriangulation
+  trian::AppendedTriangulation,
+  fe_space :: FESpace
 )
-  a = DifferentiableTriangulation(trian.a)
-  b = DifferentiableTriangulation(trian.b)
+  a = DifferentiableTriangulation(trian.a,fe_space)
+  b = DifferentiableTriangulation(trian.b,fe_space)
   return AppendedTriangulation(a,b)
 end
 
-update_trian!(trian::Triangulation,φh) = trian
+update_trian!(trian::Triangulation,U,φh) = trian
 
-function update_trian!(trian::AppendedTriangulation,φh)
-  update_trian!(trian.a,φh)
-  update_trian!(trian.b,φh)
+function update_trian!(trian::AppendedTriangulation,U,φh)
+  update_trian!(trian.a,U,φh)
+  update_trian!(trian.b,U,φh)
   return trian
 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 ? update_trian!(trian,cf) : update_trian!(trian,nothing)
+    update_trian!(trian,U,cf)
     cell_grad = f(cf)
-    update_trian!(trian,nothing) # TODO: experimental
+    update_trian!(trian,U,nothing) # TODO: experimental
     get_contribution(cell_grad,trian)
   end
   g
@@ -412,14 +411,4 @@ function FESpaces._compute_cell_ids(uh,ttrian::AppendedTriangulation)
   ids_a = FESpaces._compute_cell_ids(uh,ttrian.a)
   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
+end

src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl

@@ -11,14 +11,12 @@ Burman et al. (2017). DOI: `10.1016/j.cma.2017.09.005`.
 - `space::B`: Level-set FE space
 - `assembler::Assembler`: FE assembler
 - `params::D`: Tuple of stabilisation parameter `γg`, mesh size `h`, and
-  max steps `max_steps`
+  max steps `max_steps`, and background mesh skeleton parameters 
 - `cache`: Cache for evolver, initially `nothing`.
 
 # Note
 - The stepsize `dt = 0.1` in `RungeKutta` is a place-holder and is updated using
   the `γ` passed to `solve!`.
-- We expect the EmbeddedCollection `Ωs` to contain `:F_act` and `:dF_act`. If 
-  this is not available we add it to the recipe list in `Ωs` and a warning will appear. 
 """
 mutable struct CutFEMEvolve{A,B,C} <: Evolver
   ode_solver::ODESolver
@@ -35,56 +33,26 @@ mutable struct CutFEMEvolve{A,B,C} <: Evolver
       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 !((: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 ':Γg' under a different
-      name to avoid additional computation for cutting."
-      function F_act_recipe(cutgeo)
-        Γg = GhostSkeleton(cutgeo)
-        (; 
-          :Γ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)
+    Γg = SkeletonTriangulation(get_triangulation(dΩ_bg.quad))
+    dΓg = Measure(Γg,2get_order(V_φ))
+    n_Γg = get_normal_vector(Γg)
+    params = (;γg,h,max_steps,dΓg,n_Γg)
     new{A,B,typeof(params)}(ode_solver,Ωs,dΩ_bg,V_φ,assembler,params,nothing)
   end
 end
 
 get_ode_solver(s::CutFEMEvolve) = s.ode_solver
 get_assembler(s::CutFEMEvolve) = s.assembler
 get_space(s::CutFEMEvolve) = s.space
-get_embedded_collection(s::CutFEMEvolve) = s.Ωs
 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;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
+function get_transient_operator(φh,velh,s::CutFEMEvolve)
+  V_φ, dΩ_bg, assembler, params = s.space, s.dΩ_bg, s.assembler, s.params
+  γg, h, dΓg, n_Γg = params.γg, params.h, params.dΓg, params.n_Γg
   ϵ = 1e-20
 
-  # 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_Γg)*jump(∇(v) ⋅ n_Γg))dΓg
   mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ_bg
@@ -112,7 +80,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;use_full_skeleton=true)
+  ode_op = get_transient_operator(φh,velh,s)
   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)