Merge branch 'feature-gridap-embedded-compat' of github.com:zjwegert/…

…GridapTopOpt.jl into feature-gridap-embedded-compat

scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit.jl → scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/reinit.jl

@@ -1,69 +1,69 @@
-# Based on:
-# @article{
-#   Mallon_Thornton_Hill_Badia,
-#   title={NEURAL LEVEL SET TOPOLOGY OPTIMIZATION USING UNFITTED FINITE ELEMENTS},
-#   author={Mallon, Connor N and Thornton, Aaron W and Hill, Matthew R and Badia, Santiago},
-#   language={en}
-# }
-
-using GridapTopOpt
-
-using Gridap
-
-using GridapEmbedded
-using GridapEmbedded.LevelSetCutters
-using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData
-import Gridap.Geometry: get_node_coordinates, collect1d
-
-include("../../differentiable_trians.jl")
-
-order = 1
-n = 101
-
-_model = CartesianDiscreteModel((0,1,0,1),(n,n))
-cd = Gridap.Geometry.get_cartesian_descriptor(_model)
-h = maximum(cd.sizes)
-
-model = simplexify(_model)
-Ω = Triangulation(model)
-dΩ = Measure(Ω,2*order)
-reffe_scalar = ReferenceFE(lagrangian,Float64,order)
-V_φ = TestFESpace(model,reffe_scalar)
-
-# φh = interpolate(x->(x[1]-0.5)^2+(x[2]-0.5)^2-0.25^2,V_φ)
-φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])-0.11,V_φ)
-geo = DiscreteGeometry(φh,model)
-cutgeo = cut(model,geo)
-Γ = EmbeddedBoundary(cutgeo)
-dΓ = Measure(Γ,2*order)
-
-bgcell_to_inoutcut = compute_bgcell_to_inoutcut(model,geo);
-
-begin
-  γd = 20
-  γg = 0.1
-  ν  = 1
-  cₐ = 0.5 # <- 3 in connor's paper
-  ϵ = 1e-20
-  sgn(ϕ₀) = sign ∘ ϕ₀
-  d1(∇u) = 1 / ( ϵ + norm(∇u) )
-  W(u) =  sgn(u) * ∇(u) * (d1 ∘ (∇(u)))
-  νₐ(w) = cₐ*h * (sqrt∘( w ⋅ w ))
-  a_ν(w,u,v) = ∫((γd/h)*v*u)dΓ + ∫(νₐ(W(w))*∇(u)⋅∇(v) +  v*W(w)⋅∇(u))dΩ
-  b_ν(w,v) = ∫( sgn(w)*v )dΩ
-  res(u,v)    = a_ν(u,u,v)  - b_ν(u,v)
-  jac(u,du,v)  = a_ν(u,du,v)
-
-  op = FEOperator(res,jac,V_φ,V_φ)
-  ls = LUSolver()
-  nls = NLSolver(ftol=1e-14, iterations= 50, show_trace=true)
-  solver = FESolver(nls)
-  φh_new = FEFunction(V_φ,copy(φh.free_values))
-  Gridap.solve!(φh_new,nls,op)
-
-  writevtk(
-    Ω,"results/test_reinit",
-    cellfields=["φh"=>φh,"|∇φh|"=>norm ∘ ∇(φh),"φh_new"=>φh_new,"|∇φh_new|"=>norm ∘∇(φh_new)],
-    celldata=["inoutcut"=>bgcell_to_inoutcut]
-  )
+# Based on:
+# @article{
+#   Mallon_Thornton_Hill_Badia,
+#   title={NEURAL LEVEL SET TOPOLOGY OPTIMIZATION USING UNFITTED FINITE ELEMENTS},
+#   author={Mallon, Connor N and Thornton, Aaron W and Hill, Matthew R and Badia, Santiago},
+#   language={en}
+# }
+
+using GridapTopOpt
+
+using Gridap
+
+using GridapEmbedded
+using GridapEmbedded.LevelSetCutters
+using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData
+import Gridap.Geometry: get_node_coordinates, collect1d
+
+include("../../../differentiable_trians.jl")
+
+order = 1
+n = 101
+
+_model = CartesianDiscreteModel((0,1,0,1),(n,n))
+cd = Gridap.Geometry.get_cartesian_descriptor(_model)
+h = maximum(cd.sizes)
+
+model = simplexify(_model)
+Ω = Triangulation(model)
+dΩ = Measure(Ω,2*order)
+reffe_scalar = ReferenceFE(lagrangian,Float64,order)
+V_φ = TestFESpace(model,reffe_scalar)
+
+# φh = interpolate(x->(x[1]-0.5)^2+(x[2]-0.5)^2-0.25^2,V_φ)
+φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])-0.11,V_φ)
+geo = DiscreteGeometry(φh,model)
+cutgeo = cut(model,geo)
+Γ = EmbeddedBoundary(cutgeo)
+dΓ = Measure(Γ,2*order)
+
+bgcell_to_inoutcut = compute_bgcell_to_inoutcut(model,geo);
+
+begin
+  γd = 20
+  γg = 0.1
+  ν  = 1
+  cₐ = 0.5 # <- 3 in connor's paper
+  ϵ = 1e-20
+  sgn(ϕ₀) = sign ∘ ϕ₀
+  d1(∇u) = 1 / ( ϵ + norm(∇u) )
+  W(u) =  sgn(u) * ∇(u) * (d1 ∘ (∇(u)))
+  νₐ(w) = cₐ*h * (sqrt∘( w ⋅ w ))
+  a_ν(w,u,v) = ∫((γd/h)*v*u)dΓ + ∫(νₐ(W(w))*∇(u)⋅∇(v) +  v*W(w)⋅∇(u))dΩ
+  b_ν(w,v) = ∫( sgn(w)*v )dΩ
+  res(u,v)    = a_ν(u,u,v)  - b_ν(u,v)
+  jac(u,du,v)  = a_ν(u,du,v)
+
+  op = FEOperator(res,jac,V_φ,V_φ)
+  ls = LUSolver()
+  nls = NLSolver(ftol=1e-14, iterations= 50, show_trace=true)
+  solver = FESolver(nls)
+  φh_new = FEFunction(V_φ,copy(φh.free_values))
+  Gridap.solve!(φh_new,nls,op)
+
+  writevtk(
+    Ω,"results/test_reinit",
+    cellfields=["φh"=>φh,"|∇φh|"=>norm ∘ ∇(φh),"φh_new"=>φh_new,"|∇φh_new|"=>norm ∘∇(φh_new)],
+    celldata=["inoutcut"=>bgcell_to_inoutcut]
+  )
 end

scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unfitted_evolution.jl → scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/unfitted_evolution.jl

@@ -1,119 +1,119 @@
-struct UnfittedFEEvolution <: GridapTopOpt.LevelSetEvolution
-  spaces
-  params
-  ode_solver
-  reinit_nls
-  function UnfittedFEEvolution(Ut_φ, V_φ, dΩ, h;
-      ode_ls = LUSolver(),
-      ode_nl = NLSolver(ode_ls, show_trace=true, method=:newton, iterations=10),
-      ode_solver = RungeKutta(ode_nl, ode_ls, 0.01, :DIRK_CrankNicolson_2_2),
-      NT = 10,
-      c = 0.1,
-      reinit_nls = NLSolver(ftol=1e-14, iterations = 50, show_trace=true)
-    )
-    spaces = (Ut_φ,V_φ)
-    params = (;NT,h,c,dΩ)
-    return new(spaces,params,ode_solver,reinit_nls)
-  end
-end
-
-function GridapTopOpt.evolve!(s::UnfittedFEEvolution,φ::AbstractArray,args...)
-  _, V_φ = s.spaces
-  evolve!(s,FEFunction(V_φ,φ),args...)
-end
-
-# Based on:
-# @article{
-#   Burman_Elfverson_Hansbo_Larson_Larsson_2018,
-#   title={Shape optimization using the cut finite element method},
-#   volume={328},
-#   ISSN={00457825},
-#   DOI={10.1016/j.cma.2017.09.005},
-#   journal={Computer Methods in Applied Mechanics and Engineering},
-#   author={Burman, Erik and Elfverson, Daniel and Hansbo, Peter and Larson, Mats G. and Larsson, Karl},
-#   year={2018},
-#   month=jan,
-#   pages={242–261},
-#   language={en}
-# }
-function GridapTopOpt.evolve!(s::UnfittedFEEvolution,φh,vel,γ)
-  Ut_φ, V_φ = s.spaces
-  params = s.params
-  NT, h, c, dΩ = params.NT, params.h, params.c, params.dΩ
-  # ode_solver = s.ode_solver
-  Tf = γ*NT
-  velh = FEFunction(V_φ,vel) # <- needs to be normalised by Hilbertian projection norm
-
-  # This is temp as can't update γ in ode_solver yet
-  ode_ls = LUSolver()
-  ode_nl = NLSolver(ode_ls, show_trace=true, method=:newton, iterations=10)
-  ode_solver = RungeKutta(ode_nl, ode_ls, γ, :DIRK_CrankNicolson_2_2)
-
-  ϵ = 1e-20
-
-  # geo = DiscreteGeometry(φh,model)
-  # F = EmbeddedFacetDiscretization(LevelSetCutter(),model,geo)
-  FΓ = SkeletonTriangulation(get_triangulation(φh))#F)
-  dFΓ = Measure(FΓ,2*order)
-  n = get_normal_vector(FΓ)
-
-  d1(∇u) = 1/(ϵ + norm(∇u))
-  _n(∇u) = ∇u/(ϵ + norm(∇u))
-  β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh))
-  stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ + ∫(c*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dFΓ
-  mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ
-  forcing(t,v) = ∫(0v)dΩ
-
-  op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ,constant_forms=(true,true))
-  uht = solve(ode_solver,op,0.0,Tf,φh)
-  for (t,uh) in uht
-    if t ≈ Tf
-      copyto!(get_free_dof_values(φh),get_free_dof_values(uh))
-    end
-  end
-end
-
-function GridapTopOpt.reinit!(s::UnfittedFEEvolution,φ::AbstractArray,args...)
-  _, V_φ = s.spaces
-  reinit!(s,FEFunction(V_φ,φ),args...)
-end
-
-# Based on:
-# @article{
-#   Mallon_Thornton_Hill_Badia,
-#   title={NEURAL LEVEL SET TOPOLOGY OPTIMIZATION USING UNFITTED FINITE ELEMENTS},
-#   author={Mallon, Connor N and Thornton, Aaron W and Hill, Matthew R and Badia, Santiago},
-#   language={en}
-# }
-function GridapTopOpt.reinit!(s::UnfittedFEEvolution,φh,args...)
-  Ut_φ, V_φ = s.spaces
-  params = s.params
-  NT, h, c, dΩ = params.NT, params.h, params.c, params.dΩ
-  reinit_nls = s.reinit_nls
-
-  γd = 20
-  cₐ = 3.0 # 0.5 # <- 3 in connor's paper
-  ϵ = 1e-20
-
-  # Tmp
-  geo = DiscreteGeometry(φh,model)
-  cutgeo = cut(model,geo)
-  Γ = EmbeddedBoundary(cutgeo)
-  dΓ = Measure(Γ,2*order)
-
-  sgn(ϕ₀) = sign ∘ ϕ₀
-  d1(∇u) = 1 / ( ϵ + norm(∇u) )
-  W(u) =  sgn(u) * ∇(u) * (d1 ∘ (∇(u)))
-  νₐ(w) = cₐ*h*(sqrt∘( w ⋅ w ))
-  a_ν(w,u,v) = ∫((γd/h)*v*u)dΓ + ∫(νₐ(W(w))*∇(u)⋅∇(v) + v*W(w)⋅∇(u))dΩ
-  b_ν(w,v) = ∫( sgn(w)*v )dΩ
-  res(u,v)    = a_ν(u,u,v)  - b_ν(u,v)
-  jac(u,du,v)  = a_ν(u,du,v)
-
-  op = FEOperator(res,jac,V_φ,V_φ)
-  Gridap.solve!(φh,reinit_nls,op)
-end
-
-function GridapTopOpt.get_dof_Δ(s::UnfittedFEEvolution)
-  s.params.h
+struct UnfittedFEEvolution <: GridapTopOpt.LevelSetEvolution
+  spaces
+  params
+  ode_solver
+  reinit_nls
+  function UnfittedFEEvolution(Ut_φ, V_φ, dΩ, h;
+      ode_ls = LUSolver(),
+      ode_nl = NLSolver(ode_ls, show_trace=true, method=:newton, iterations=10),
+      ode_solver = RungeKutta(ode_nl, ode_ls, 0.01, :DIRK_CrankNicolson_2_2),
+      NT = 10,
+      c = 0.1,
+      reinit_nls = NLSolver(ftol=1e-14, iterations = 50, show_trace=true)
+    )
+    spaces = (Ut_φ,V_φ)
+    params = (;NT,h,c,dΩ)
+    return new(spaces,params,ode_solver,reinit_nls)
+  end
+end
+
+function GridapTopOpt.evolve!(s::UnfittedFEEvolution,φ::AbstractArray,args...)
+  _, V_φ = s.spaces
+  evolve!(s,FEFunction(V_φ,φ),args...)
+end
+
+# Based on:
+# @article{
+#   Burman_Elfverson_Hansbo_Larson_Larsson_2018,
+#   title={Shape optimization using the cut finite element method},
+#   volume={328},
+#   ISSN={00457825},
+#   DOI={10.1016/j.cma.2017.09.005},
+#   journal={Computer Methods in Applied Mechanics and Engineering},
+#   author={Burman, Erik and Elfverson, Daniel and Hansbo, Peter and Larson, Mats G. and Larsson, Karl},
+#   year={2018},
+#   month=jan,
+#   pages={242–261},
+#   language={en}
+# }
+function GridapTopOpt.evolve!(s::UnfittedFEEvolution,φh,vel,γ)
+  Ut_φ, V_φ = s.spaces
+  params = s.params
+  NT, h, c, dΩ = params.NT, params.h, params.c, params.dΩ
+  # ode_solver = s.ode_solver
+  Tf = γ*NT
+  velh = FEFunction(V_φ,vel) # <- needs to be normalised by Hilbertian projection norm
+
+  # This is temp as can't update γ in ode_solver yet
+  ode_ls = LUSolver()
+  ode_nl = NLSolver(ode_ls, show_trace=true, method=:newton, iterations=10)
+  ode_solver = RungeKutta(ode_nl, ode_ls, γ, :DIRK_CrankNicolson_2_2)
+
+  ϵ = 1e-20
+
+  # geo = DiscreteGeometry(φh,model)
+  # F = EmbeddedFacetDiscretization(LevelSetCutter(),model,geo)
+  FΓ = SkeletonTriangulation(get_triangulation(φh))#F)
+  dFΓ = Measure(FΓ,2*order)
+  n = get_normal_vector(FΓ)
+
+  d1(∇u) = 1/(ϵ + norm(∇u))
+  _n(∇u) = ∇u/(ϵ + norm(∇u))
+  β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh))
+  stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ + ∫(c*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dFΓ
+  mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ
+  forcing(t,v) = ∫(0v)dΩ
+
+  op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ,constant_forms=(true,true))
+  uht = solve(ode_solver,op,0.0,Tf,φh)
+  for (t,uh) in uht
+    if t ≈ Tf
+      copyto!(get_free_dof_values(φh),get_free_dof_values(uh))
+    end
+  end
+end
+
+function GridapTopOpt.reinit!(s::UnfittedFEEvolution,φ::AbstractArray,args...)
+  _, V_φ = s.spaces
+  reinit!(s,FEFunction(V_φ,φ),args...)
+end
+
+# Based on:
+# @article{
+#   Mallon_Thornton_Hill_Badia,
+#   title={NEURAL LEVEL SET TOPOLOGY OPTIMIZATION USING UNFITTED FINITE ELEMENTS},
+#   author={Mallon, Connor N and Thornton, Aaron W and Hill, Matthew R and Badia, Santiago},
+#   language={en}
+# }
+function GridapTopOpt.reinit!(s::UnfittedFEEvolution,φh,args...)
+  Ut_φ, V_φ = s.spaces
+  params = s.params
+  NT, h, c, dΩ = params.NT, params.h, params.c, params.dΩ
+  reinit_nls = s.reinit_nls
+
+  γd = 20
+  cₐ = 3.0 # 0.5 # <- 3 in connor's paper
+  ϵ = 1e-20
+
+  # Tmp
+  geo = DiscreteGeometry(φh,model)
+  cutgeo = cut(model,geo)
+  Γ = EmbeddedBoundary(cutgeo)
+  dΓ = Measure(Γ,2*order)
+
+  sgn(ϕ₀) = sign ∘ ϕ₀
+  d1(∇u) = 1 / ( ϵ + norm(∇u) )
+  W(u) =  sgn(u) * ∇(u) * (d1 ∘ (∇(u)))
+  νₐ(w) = cₐ*h*(sqrt∘( w ⋅ w ))
+  a_ν(w,u,v) = ∫((γd/h)*v*u)dΓ + ∫(νₐ(W(w))*∇(u)⋅∇(v) + v*W(w)⋅∇(u))dΩ
+  b_ν(w,v) = ∫( sgn(w)*v )dΩ
+  res(u,v)    = a_ν(u,u,v)  - b_ν(u,v)
+  jac(u,du,v)  = a_ν(u,du,v)
+
+  op = FEOperator(res,jac,V_φ,V_φ)
+  Gridap.solve!(φh,reinit_nls,op)
+end
+
+function GridapTopOpt.get_dof_Δ(s::UnfittedFEEvolution)
+  s.params.h
 end