added unfitted evolve and reinit test scripts

scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl

@@ -0,0 +1,83 @@
+# 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}
+# }
+
+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->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) # <- Already SDF
+# φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])-0.11,V_φ) # This needs to reinitialised, can use `reinit_based_on_connor_paper.jl`
+
+geo = DiscreteGeometry(φh,model)
+cutgeo = cut(model,geo)
+Γ = EmbeddedBoundary(cutgeo)
+dΓ = Measure(Γ,2*order)
+
+F = EmbeddedFacetDiscretization(LevelSetCutter(),model,geo)
+FΓ = SkeletonTriangulation(F)
+dFΓ = Measure(FΓ,2*order)
+
+dt = 0.01
+Tn = 5*dt
+sysslvr_l = LUSolver()
+sysslvr_nl = NLSolver(sysslvr_l, show_trace=true, method=:newton, iterations=10)
+ode_solver = RungeKutta(sysslvr_nl, sysslvr_l, dt, :DIRK_CrankNicolson_2_2)
+
+c = 0.1
+
+Ut_φ = TransientTrialFESpace(V_φ,t -> (x->-1))
+n = get_normal_vector(FΓ)
+velh = interpolate(x->-1,V_φ)
+d1(∇u) = 1 / ( ϵ + norm(∇u) )
+_n(∇u) = ∇u/(10^-20+norm(∇u))
+β = velh*∇(φh)/(10^-20+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,Tn,φh)
+for (t,uh) in uht
+  if t ≈ Tn
+    copyto!(get_free_dof_values(φh),get_free_dof_values(uh))
+  end
+end
+
+# You can compare the zero-iso curve generated by the following
+φh_expected_contour = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25+Tn,V_φ)
+
+writevtk(
+  Ω,"results/test_evolve",
+  cellfields=["φh"=>φh,"φh_expected_contour"=>φh_expected_contour,"|∇φh|"=>norm ∘ ∇(φh)]
+)

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

@@ -1,3 +1,11 @@
+# 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
@@ -7,7 +15,7 @@ using GridapEmbedded.LevelSetCutters
 using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData
 import Gridap.Geometry: get_node_coordinates, collect1d
 
-include("../differentiable_trians.jl")
+include("../../differentiable_trians.jl")
 
 order = 1
 n = 101

scripts/Embedded/MWEs/zach_embedded.jl

@@ -14,8 +14,8 @@ n = 101
 N = 16
 a = 1
 begin
-model = CartesianDiscreteModel((0,1,0,1),(n,n))
-model = simplexify(model)
+_model = CartesianDiscreteModel((0,1,0,1),(n,n))
+model = simplexify(_model)
 Ω = Triangulation(model)
 dΩ = Measure(Ω,2*order)
 
@@ -28,8 +28,8 @@ V_φ = TestFESpace(model,reffe_scalar)
 # φh = interpolate(x->(1+0.25)abs(x[1]-0.5)+0abs(x[2]-0.5)-0.25,V_φ) # Straight lines with scaling
 # φh = interpolate(x->abs(x[1]-0.5)+0abs(x[2]-0.5)-0.25/(1+0.25),V_φ) # Straight lines without scaling
 # φh = interpolate(x->tan(-pi/3)*(x[1]-0.5)+(x[2]-0.5),V_φ) # Angled interface
-φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.25,V_φ) # Circle
-# φh = interpolate(x->a*(cos(2π*x[1])*cos(2π*x[2])-0.11),V_φ) # "Regular" LSF
+# φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.25,V_φ) # Circle
+φh = interpolate(x->a*(cos(2π*x[1])*cos(2π*x[2])-0.11),V_φ) # "Regular" LSF
 x_φ = get_free_dof_values(φh)
 
 if any(isapprox(0.0;atol=10^-10),x_φ)
@@ -52,20 +52,20 @@ diff_Ωout = DifferentiableTriangulation(Ωout)
 fh = interpolate(x->x[1]+x[2],V_φ)
 
 dΩin = Measure(diff_Ωin,2*order)
-j_in(φ) = ∫(cos ∘ fh)dΩin
+j_in(φ) = ∫(fh)dΩin
 dj_in = gradient(j_in,φh)
 dj_vec_in = assemble_vector(dj_in,V_φ)
 norm(dj_vec_in)
 
 dΩout = Measure(diff_Ωout,2*order)
-j_out(φ) = ∫(cos ∘ fh)dΩout
+j_out(φ) = ∫(fh)dΩout
 dj_out = gradient(j_out,φh)
 dj_vec_out = -assemble_vector(dj_out,V_φ)
 norm(dj_vec_out)
 
 Γ = EmbeddedBoundary(cutgeo)
 dΓ = Measure(Γ,2*order)
-dj_expected(q) = ∫(-(cos ∘ fh)*q/(norm ∘ (∇(φh))))dΓ
+dj_expected(q) = ∫(-(fh)*q/(norm ∘ (∇(φh))))dΓ
 dj_exp_vec = assemble_vector(dj_expected,V_φ)
 norm(dj_exp_vec)
 
@@ -108,6 +108,7 @@ V_φ = TestFESpace(model,reffe_scalar)
 # φh = interpolate(x->tan(-pi/3)*(x[1]-0.5)+(x[2]-0.5),V_φ) # Angled interface
 # φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2+(x[3]-0.5)^2)-0.25,V_φ) # Sphere
 φh = interpolate(x->a*(cos(2π*x[1])*cos(2π*x[2])*cos(2π*x[3])-0.11),V_φ) # "Regular" LSF
+# φh = interpolate(x->0.3(cos(2π*(x[1]-0.5))*cos(2π*(x[2]-0.5))*cos(2π*(x[3]-0.5)))+0.2(cos(2π*(x[1]-0.5))+cos(2π*(x[2]-0.5))+cos(2π*(x[3]-0.5)))+0.1(cos(2π*2*(x[1]-0.5))*cos(2π*2*(x[2]-0.5))*cos(2π*2*(x[3]-0.5)))+0.1(cos(2π*2*(x[1]-0.5))+cos(2π*2*(x[2]-0.5))+cos(2π*2*(x[3]-0.5)))+0.05(cos(2π*3*(x[1]-0.5))+cos(2π*3*(x[2]-0.5))+cos(2π*3*(x[3]-0.5)))+0.1(cos(2π*(x[1]-0.5))*cos(2π*(x[2]-0.5))+cos(2π*(x[2]-0.5))*cos(2π*(x[3]-0.5))+cos(2π*(x[3]-0.5))*cos(2π*(x[1]-0.5))),V_φ)
 x_φ = get_free_dof_values(φh)
 
 if any(isapprox(0.0;atol=10^-10),x_φ)
@@ -130,7 +131,7 @@ diff_Ωout = DifferentiableTriangulation(Ωout)
 fh = interpolate(x->x[1]+x[2]+x[3],V_φ)
 
 dΩin = Measure(diff_Ωin,2*order)
-j_in(φ) = ∫(fh)dΩin # <- Cannot take fh due to missing method?
+j_in(φ) = ∫(fh)dΩin
 dj_in = gradient(j_in,φh)
 dj_vec_in = assemble_vector(dj_in,V_φ)
 norm(dj_vec_in)