Additional test

scripts/Embedded/2d_elastic_compliance.jl

@@ -0,0 +1,107 @@
+using Pkg; Pkg.activate()
+
+using Gridap,GridapTopOpt
+include("embedded_measures.jl")
+
+function main(path="./results/UnfittedFEM_elastic_compliance_ALM/")
+  ## Parameters
+  order = 1
+  xmax,ymax=(2.0,1.0)
+  prop_Γ_N = 0.2
+  dom = (0,xmax,0,ymax)
+  el_size = (200,100)
+  γ = 0.1
+  γ_reinit = 0.5
+  max_steps = floor(Int,order*minimum(el_size)/10)
+  tol = 1/(5order^2)/minimum(el_size)
+  C = isotropic_elast_tensor(2,1.,0.3)
+  η_coeff = 2
+  α_coeff = 4max_steps*γ
+  vf = 0.4
+  g = VectorValue(0,-1)
+  iter_mod = 1
+  mkpath(path)
+
+  ## FE Setup
+  model = CartesianDiscreteModel(dom,el_size)
+  el_Δ = get_el_Δ(model)
+  f_Γ_D(x) = (x[1] ≈ 0.0)
+  f_Γ_N(x) = (x[1] ≈ xmax && ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/2 + 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Ω)
+
+  ## Spaces
+  reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
+  reffe_scalar = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"])
+  U = TrialFESpace(V,VectorValue(0.0,0.0))
+  V_φ = TestFESpace(model,reffe_scalar)
+  V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"])
+  U_reg = TrialFESpace(V_reg,0)
+
+  ## Create FE functions
+  φh = interpolate(initial_lsf(4,0.2),V_φ)
+  embedded_meas = EmbeddedMeasureCache(φh,V_φ)
+  update_meas(φ) = update_embedded_measures!(φ,embedded_meas)
+  get_meas(φ) = get_embedded_measures(φ,embedded_meas)
+
+  ## Interpolation and weak form
+  interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
+  I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
+
+  a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(C ⊙ ε(u) ⊙ ε(v))dΩ1 + ∫((10^-6*C) ⊙ ε(u) ⊙ ε(v))dΩ2
+  l(v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(v⋅g)dΓ_N
+
+  ## Optimisation functionals
+  J(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(C ⊙ ε(u) ⊙ ε(u))dΩ1
+  dJ(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫((C ⊙ ε(u) ⊙ ε(u))*q)dΓ;
+  Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf/vol_D)dΩ;
+  dVol(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(-1/vol_D*q)dΓ
+
+  ## IntegrandWithEmbeddedMeasure
+  a_iem = IntegrandWithEmbeddedMeasure(a,(dΩ,dΓ_N),update_meas)
+  l_iem = IntegrandWithEmbeddedMeasure(l,(dΩ,dΓ_N),get_meas)
+  J_iem = IntegrandWithEmbeddedMeasure(J,(dΩ,dΓ_N),get_meas)
+  dJ_iem = IntegrandWithEmbeddedMeasure(dJ,(dΩ,dΓ_N),get_meas)
+  Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),get_meas)
+  dVol_iem = IntegrandWithEmbeddedMeasure(dVol,(dΩ,dΓ_N),get_meas)
+
+  ## Finite difference solver and level set function
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
+
+  ## Setup solver and FE operators
+  state_map = AffineFEStateMap(a_iem,l_iem,U,V,V_φ,U_reg,φh,(dΩ,dΓ_N))
+  pcfs = PDEConstrainedFunctionals(J_iem,[Vol_iem],state_map,analytic_dJ=dJ_iem,analytic_dC=[dVol_iem])
+
+  ## 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)
+
+  ## Optimiser
+  rm(path,force=true,recursive=true)
+  mkpath(path)
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;reinit_mod=5,
+    γ,γ_reinit,verbose=true,constraint_names=[:Vol])
+  for (it,uh,φh) in optimiser
+    _Ω1,_,_Γ = get_embedded_triangulations(embedded_meas)
+    if iszero(it % iter_mod)
+      writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+      writevtk(_Γ,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)
+  _Ω1,_,_Γ = get_embedded_triangulations(embedded_meas)
+  writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+  writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)])
+end
+
+main()