From 38190c24da77c8b26c11152565552710f22d630d Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Mon, 16 Sep 2024 20:05:20 +1000 Subject: [PATCH] added unfitted evolve and reinit test scripts --- .../MWEs/Unfitted_Evolve&Reinit/evolve.jl | 83 +++++++++++++++++++ .../reinit.jl} | 10 ++- scripts/Embedded/MWEs/zach_embedded.jl | 17 ++-- 3 files changed, 101 insertions(+), 9 deletions(-) create mode 100644 scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl rename scripts/Embedded/MWEs/{reinit_based_on_connor_paper.jl => Unfitted_Evolve&Reinit/reinit.jl} (82%) diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl new file mode 100644 index 00000000..8d5ced03 --- /dev/null +++ b/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)] +) \ No newline at end of file diff --git a/scripts/Embedded/MWEs/reinit_based_on_connor_paper.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit.jl similarity index 82% rename from scripts/Embedded/MWEs/reinit_based_on_connor_paper.jl rename to scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit.jl index 0ef58568..f12f9afe 100644 --- a/scripts/Embedded/MWEs/reinit_based_on_connor_paper.jl +++ b/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 diff --git a/scripts/Embedded/MWEs/zach_embedded.jl b/scripts/Embedded/MWEs/zach_embedded.jl index 94e899ee..9d5a049c 100644 --- a/scripts/Embedded/MWEs/zach_embedded.jl +++ b/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)