diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/reinit.jl similarity index 94% rename from scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit.jl rename to scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/reinit.jl index f12f9afe..b3b9b4a1 100644 --- a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit.jl +++ b/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 \ No newline at end of file diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unfitted_evolution.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/unfitted_evolution.jl similarity index 96% rename from scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unfitted_evolution.jl rename to scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/unfitted_evolution.jl index 20aac60f..9d8c5f60 100644 --- a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unfitted_evolution.jl +++ b/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 \ No newline at end of file diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unitted_evolution_test.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/unitted_evolution_test.jl similarity index 95% rename from scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unitted_evolution_test.jl rename to scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/unitted_evolution_test.jl index 56460f84..3184e2cb 100644 --- a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unitted_evolution_test.jl +++ b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/unitted_evolution_test.jl @@ -1,66 +1,66 @@ -using GridapTopOpt - -using Gridap - -using GridapEmbedded -using GridapEmbedded.LevelSetCutters -using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.ODEs, Gridap.Adaptivity -import Gridap.Geometry: get_node_coordinates, collect1d - -include("unfitted_evolution.jl") -include("../../differentiable_trians.jl") - -order = 1 -n = 101 - -_model = CartesianDiscreteModel((0,1,0,1),(n,n)) -cd = Gridap.Geometry.get_cartesian_descriptor(_model) -base_model = UnstructuredDiscreteModel(_model) -ref_model = refine(base_model, refinement_method = "barycentric") -model = ref_model.model -# model = simplexify(_model) -h = maximum(cd.sizes) - -Ω = Triangulation(model) -dΩ = Measure(Ω,2*order) -reffe_scalar = ReferenceFE(lagrangian,Float64,order) -V_φ = TestFESpace(model,reffe_scalar) -Ut_φ = TransientTrialFESpace(V_φ) - -α = 2.0*h -a_hilb(p,q) = ∫(α^2*∇(p)⋅∇(q) + p*q)dΩ; -vel_ext = VelocityExtension(a_hilb,V_φ,V_φ) - -# φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) # <- Already SDF -φh = interpolate(x->cos(4π*x[1])*cos(4π*x[2])-0.11,V_φ) - -function compute_test_velh(φh) - geo = DiscreteGeometry(φh,model) - cutgeo = cut(model,geo) - Ωin = Triangulation(cutgeo,PHYSICAL_IN).a - diff_Ωin = DifferentiableTriangulation(Ωin) - dΩin = Measure(diff_Ωin,2*order) - j_in(φ) = ∫(φ)dΩin - dj_in = gradient(j_in,φh) - dj_vec_in = assemble_vector(dj_in,V_φ) - project!(vel_ext,dj_vec_in) - vel = dj_vec_in/sqrt(dj_vec_in'*vel_ext.K*dj_vec_in) - return FEFunction(V_φ,vel) -end - -ls_evo = UnfittedFEEvolution(Ut_φ, V_φ, dΩ, h; NT = 5, c = 0.1, - reinit_nls = NLSolver(ftol=1e-12, iterations = 50, show_trace=true)) -reinit!(ls_evo,φh) # <- inital sdf - -velh = compute_test_velh(φh) -# velh = interpolate(x->-1,V_φ) -evolve!(ls_evo,φh,get_free_dof_values(velh),0.01) -reinit!(ls_evo,φh) - -# You can compare the zero-iso curve for advected circle via -# φh_expected_contour = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25+ls_evo.params.NT*0.01,V_φ) - -writevtk( - Ω,"results/test_evolve", - cellfields=["φh"=>φh,"|∇φh|"=>norm ∘ ∇(φh),"velh"=>velh]#,"φh_expected_contour"=>φh_expected_contour] +using GridapTopOpt + +using Gridap + +using GridapEmbedded +using GridapEmbedded.LevelSetCutters +using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.ODEs, Gridap.Adaptivity +import Gridap.Geometry: get_node_coordinates, collect1d + +include("unfitted_evolution.jl") +include("../../../differentiable_trians.jl") + +order = 1 +n = 101 + +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +cd = Gridap.Geometry.get_cartesian_descriptor(_model) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model +# model = simplexify(_model) +h = maximum(cd.sizes) + +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order) +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_φ = TestFESpace(model,reffe_scalar) +Ut_φ = TransientTrialFESpace(V_φ) + +α = 2.0*h +a_hilb(p,q) = ∫(α^2*∇(p)⋅∇(q) + p*q)dΩ; +vel_ext = VelocityExtension(a_hilb,V_φ,V_φ) + +# φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) # <- Already SDF +φh = interpolate(x->cos(4π*x[1])*cos(4π*x[2])-0.11,V_φ) + +function compute_test_velh(φh) + geo = DiscreteGeometry(φh,model) + cutgeo = cut(model,geo) + Ωin = Triangulation(cutgeo,PHYSICAL_IN).a + diff_Ωin = DifferentiableTriangulation(Ωin) + dΩin = Measure(diff_Ωin,2*order) + j_in(φ) = ∫(φ)dΩin + dj_in = gradient(j_in,φh) + dj_vec_in = assemble_vector(dj_in,V_φ) + project!(vel_ext,dj_vec_in) + vel = dj_vec_in/sqrt(dj_vec_in'*vel_ext.K*dj_vec_in) + return FEFunction(V_φ,vel) +end + +ls_evo = UnfittedFEEvolution(Ut_φ, V_φ, dΩ, h; NT = 5, c = 0.1, + reinit_nls = NLSolver(ftol=1e-12, iterations = 50, show_trace=true)) +reinit!(ls_evo,φh) # <- inital sdf + +velh = compute_test_velh(φh) +# velh = interpolate(x->-1,V_φ) +evolve!(ls_evo,φh,get_free_dof_values(velh),0.01) +reinit!(ls_evo,φh) + +# You can compare the zero-iso curve for advected circle via +# φh_expected_contour = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25+ls_evo.params.NT*0.01,V_φ) + +writevtk( + Ω,"results/test_evolve", + cellfields=["φh"=>φh,"|∇φh|"=>norm ∘ ∇(φh),"velh"=>velh]#,"φh_expected_contour"=>φh_expected_contour] ) \ No newline at end of file diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolution_test.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolution_test.jl new file mode 100644 index 00000000..7077ea4d --- /dev/null +++ b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolution_test.jl @@ -0,0 +1,88 @@ +using GridapTopOpt +using GridapSolvers.NonlinearSolvers + +using Gridap +using Gridap.Adaptivity + +using GridapEmbedded +using GridapEmbedded.LevelSetCutters +using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData +import Gridap.Geometry: get_node_coordinates, collect1d + +include("../../differentiable_trians.jl") + +using GridapTopOpt: LevelSetEvolution, instantiate_caches +using Gridap.Algebra: NonlinearSolver +using Gridap.ODEs: ODESolver +using GridapDistributed: DistributedFESpace +import Gridap: solve! + +include("../../../../src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl") + +order = 1 +n = 101 +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +cd = Gridap.Geometry.get_cartesian_descriptor(_model) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model +# model = simplexify(_model) +h = maximum(cd.sizes) + +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order) +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_φ = TestFESpace(model,reffe_scalar) + +α = 2.0*h +a_hilb(p,q) = ∫(α^2*∇(p)⋅∇(q) + p*q)dΩ; +vel_ext = VelocityExtension(a_hilb,V_φ,V_φ) + +φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) # <- Already SDF +# φh = interpolate(x->-(x[1]-0.5)^2-(x[2]-0.5)^2+0.25^2,V_φ) +# φh = interpolate(x->cos(4π*x[1])*cos(4π*x[2])-0.11,V_φ) + +function compute_test_velh(φh) + geo = DiscreteGeometry(φh,model) + cutgeo = cut(model,geo) + Ωin = Triangulation(cutgeo,PHYSICAL_IN).a + diff_Ωin = DifferentiableTriangulation(Ωin) + dΩin = Measure(diff_Ωin,2*order) + j_in(φ) = ∫(1)dΩin + dj_in = gradient(j_in,φh) + dj_vec_in = assemble_vector(dj_in,V_φ) + project!(vel_ext,dj_vec_in) + vel = dj_vec_in/sqrt(dj_vec_in'*vel_ext.K*dj_vec_in) + return FEFunction(V_φ,vel) +end + +ls_evo = CutFEMEvolve(model,V_φ,dΩ,h) +ls_reinit = StabilisedReinit(model,V_φ,dΩ,h) +evo = UnfittedFEEvolution(ls_evo,ls_reinit) +# reinit!(evo,φh); + +φ0 = copy(get_free_dof_values(φh)) +φh0 = FEFunction(V_φ,φ0) + +# velh = compute_test_velh(φh) +velh = interpolate(x->-1,V_φ) +evolve!(evo,φh,velh,0.1) + +ls_evo = CutFEMEvolve(model,V_φ,dΩ,h) +ls_reinit = StabilisedReinit(model,V_φ,dΩ,h) +evo = UnfittedFEEvolution(ls_evo,ls_reinit) + +velh = interpolate(x->1,V_φ) +evolve!(evo,φh,velh,0.1) + +Δt = _compute_Δt(h,0.1,get_free_dof_values(velh)) +φh_expected_contour = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25+evo.evolver.params.max_steps*Δt,V_φ) + +writevtk( + Ω,"results/test_evolve", + cellfields=[ + "φh0"=>φh0,"|∇φh0|"=>norm ∘ ∇(φh0), + "φh_advect"=>φh,"|∇φh_advect|"=>norm ∘ ∇(φh), + # "φh_expected_contour"=>φh_expected_contour + ] +) \ No newline at end of file diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl index 8d5ced03..3ac907fb 100644 --- a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl +++ b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl @@ -80,4 +80,29 @@ end writevtk( Ω,"results/test_evolve", cellfields=["φh"=>φh,"φh_expected_contour"=>φh_expected_contour,"|∇φh|"=>norm ∘ ∇(φh)] -) \ No newline at end of file +) + + +sol,state = Base.iterate(uht) +_tF = state[2][1] +_stateF = state[2][2] +_state0 = state[2][3] +_uF = state[2][4] +_odecache = state[2][5] + +_odeslvrcache, _odeopcache = _odecache + +_odeslvrcache +_odeopcache.const_forms # contains linear systems + +# For each call of evolve need to: +# 1. set _tF to zero +# 2. copy φh to state0 +# 3. reassemble matrices in _odeopcache.const_forms +# 4. solve +# 5. copy stateF to φh + +# User struct: `UnfittedFEEvolution <: LevelSetEvolution` with methods evolve! and reinit! +# Internal structs: `CutFEMEvolve <: Evolver`, `PicardReinit <: Reinitialiser` +# Types that inherit from `UnfittedFEMethod` iterate to re-use cache etc. + diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit_test.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit_test.jl new file mode 100644 index 00000000..94a077ec --- /dev/null +++ b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit_test.jl @@ -0,0 +1,56 @@ +using GridapTopOpt +using GridapSolvers.NonlinearSolvers + +using Gridap +using Gridap.Adaptivity + +using GridapEmbedded +using GridapEmbedded.LevelSetCutters +using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData +import Gridap.Geometry: get_node_coordinates, collect1d + +include("../../differentiable_trians.jl") + +using GridapTopOpt: LevelSetEvolution, instantiate_caches +using Gridap.Algebra: NonlinearSolver +using Gridap.ODEs: ODESolver +using GridapDistributed: DistributedFESpace +import Gridap: solve! + +include("../../../../src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl") + +order = 1 +n = 101 +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +cd = Gridap.Geometry.get_cartesian_descriptor(_model) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model +h = maximum(cd.sizes) + +Ω = 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_φ) + +ls_evo = CutFEMEvolve(model,V_φ,dΩ,h) +ls_reinit = StabilisedReinit(model,V_φ,dΩ,h) + +φ0 = copy(get_free_dof_values(φh)) +φh0 = FEFunction(V_φ,φ0) + +evo = UnfittedFEEvolution(ls_evo,ls_reinit) +reinit!(evo,φh); + +L2error(u) = ∫(u*u)dΩ +@assert sum(L2error((norm ∘ ∇(φh_reinit)))) - 1 < 10^-4 + +writevtk( + Ω,"results/test_evolve", + cellfields=[ + "φh0"=>φh0,"|∇φh0|"=>norm ∘ ∇(φh0), + "φh_reinit"=>φh_reinit,"|∇φh_reinit|"=>norm ∘ ∇(φh_reinit) + ] +) \ No newline at end of file diff --git a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl new file mode 100644 index 00000000..626329a9 --- /dev/null +++ b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl @@ -0,0 +1,132 @@ +""" + mutable struct CutFEMEvolve{V,M} <: Evolver + +CutFEM method for level-set evolution based on method developed by +Burman et al. (2017). DOI: `10.1016/j.cma.2017.09.005`. + +# Parameters +- `ode_solver::ODESolver`: ODE solver +- `model::A`: FE model +- `space::B`: Level-set FE space +- `dΩ::C`: Measure for integration +- `assembler::Assembler`: FE assembler +- `params::D`: Tuple of stabilisation parameter `Γg`, mesh size `h`, + max steps `max_steps`, and FE space `order` +- `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!`. +""" +mutable struct CutFEMEvolve{A,B,C,D} <: Evolver + ode_solver::ODESolver + model::A + space::B + dΩ::C + assembler::Assembler + params::D + cache + function CutFEMEvolve(model::A,V_φ::B,dΩ::C,h::Real; + max_steps=10, + Γg = 0.1, + ode_ls = LUSolver(), + 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,C} + params = (;Γg,h,max_steps,order=get_order(V_φ)) + new{A,B,C,typeof(params)}(ode_solver,model,V_φ,dΩ,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_model(s::CutFEMEvolve) = s.model +get_measure(s::CutFEMEvolve) = s.dΩ +get_params(s::CutFEMEvolve) = s.params +get_cache(s::CutFEMEvolve) = s.cache + +function get_stiffness_matrix_map(φh,velh,dΩ,model,order,Γg,h) + geo = DiscreteGeometry(φh,model) + cutgeo = cut(model,geo) + Ω_act = Triangulation(cutgeo,ACTIVE) + F_act = SkeletonTriangulation(Ω_act) + dF_act = Measure(F_act,2*order) + n = get_normal_vector(F_act) + ϵ = 1e-20 + + β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh)) + return (t,u,v) -> ∫((β ⋅ ∇(u)) * v)dΩ + ∫(Γg*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dF_act +end + +function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing) + ode_solver, model, V_φ, dΩ, assembler = s.ode_solver, s.model, s.space, s.dΩ, s.assembler + Γg, h, max_steps, order = s.params + + # Weak form + stiffness = get_stiffness_matrix_map(φh,velh,dΩ,model,order,Γg,h) + mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ + forcing(t,v) = ∫(0v)dΩ + + # Setup FE operator and solver + Ut_φ = TransientTrialFESpace(V_φ) + ode_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ; + constant_forms=(true,true),assembler) + 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) + + # March + march = Base.iterate(ode_sol) + data, state = march + while march !== nothing + data, state = march + march = Base.iterate(ode_sol,state) + end + + # Update φh and cache + _, φhF = data + copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) + s.cache = (ode_op, state) + + return φh +end + +function solve!(s::CutFEMEvolve,φh,velh,γ,cache) + ode_solver, model, dΩ, V_φ, assembler, = s.ode_solver, s.model, s.dΩ, s.space, s.assembler + Γg, h, max_steps, order = s.params + ode_op, state0 = cache + + U, (tF, stateF, state0, uF, odecache) = state0 + odeslvrcache, odeopcache = odecache + K, M = odeopcache.const_forms + + # Update odeopcache and solver + stiffness = get_stiffness_matrix_map(φh,velh,dΩ,model,order,Γg,h) + assemble_matrix!((u,v)->stiffness(0.0,u,v),K,assembler,V_φ,V_φ) + state = U, (0.0, stateF, state0, uF, odecache) + dt = _compute_Δt(h,γ,get_free_dof_values(velh)) + ode_solver.dt = dt + + # March + ode_sol = solve(ode_solver,ode_op,0.0,dt*max_steps,φh) + march = Base.iterate(ode_sol,state) + data, state = march + while march !== nothing + data, state = march + march = Base.iterate(ode_sol,state) + end + + # Update φh and cache + _, φhF = data + copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) + s.cache = (ode_op, state) + + return φh +end + +function _compute_Δt(h,γ,vel) + T = eltype(γ) + v_norm = maximum(abs,vel) + return γ * h / (eps(T)^2 + v_norm) +end \ No newline at end of file diff --git a/src/LevelSetEvolution/UnfittedEvolution/MutableRungeKutta.jl b/src/LevelSetEvolution/UnfittedEvolution/MutableRungeKutta.jl new file mode 100644 index 00000000..b9597c5e --- /dev/null +++ b/src/LevelSetEvolution/UnfittedEvolution/MutableRungeKutta.jl @@ -0,0 +1,301 @@ +import Gridap.ODEs: allocate_odecache,ode_march!,allocate_odecache, + ode_march!,_update_dimrk!,AbstractTableau,DiagonallyImplicitTableau, + AbstractQuasilinearODE,AbstractLinearODE,ODEOperator,TableauType,ExplicitTableau, + allocate_odeopcache,ODEOperatorType,is_form_constant,allocate_jacobian,allocate_residual, + NumericalSetup,update_odeopcache!,LinearStageOperator,_setindex_all!,axpy! + +######################## +# MutableDIMRungeKutta # +######################## +""" + struct MutableDIMRungeKutta <: ODESolver end + +Diagonally-implicit Runge-Kutta ODE solver. +``` +residual(tx, ux, vx) = 0, + +tx = t_n + c[i] * dt +ux = u_n + dt * ∑_{1 ≤ j < i} A[i, j] * slopes[j] + dt * A[i, i] * x +vx = x, + +u_(n+1) = u_n + dt * ∑_{1 ≤ i ≤ s} b[i] * slopes[i]. +``` +""" +mutable struct MutableDIMRungeKutta <: ODESolver + sysslvr_nl::NonlinearSolver + sysslvr_l::NonlinearSolver + dt::Real + tableau::AbstractTableau{DiagonallyImplicitTableau} +end + +################## +# Nonlinear case # +################## +function allocate_odecache( + odeslvr::MutableDIMRungeKutta, odeop::ODEOperator, + t0::Real, us0::NTuple{1,AbstractVector} +) + u0 = us0[1] + us0N = (u0, u0) + odeopcache = allocate_odeopcache(odeop, t0, us0N) + + ui_pre, ui = zero(u0), zero(u0) + num_stages = length(get_nodes(odeslvr.tableau)) + slopes = [zero(u0) for _ in 1:num_stages] + + odeoptype = ODEOperatorType(odeop) + has_explicit = odeoptype <: AbstractQuasilinearODE + is_semilinear = odeoptype <: AbstractSemilinearODE + mass_constant = is_form_constant(odeop, 1) + reuse = (is_semilinear && mass_constant) + + J, r = nothing, nothing + if has_explicit + # Allocate J, r if there are explicit stages + A = get_matrix(odeslvr.tableau) + if any(i -> iszero(A[i, i]), axes(A, 2)) + J = allocate_jacobian(odeop, t0, us0N, odeopcache) + r = allocate_residual(odeop, t0, us0N, odeopcache) + end + end + + sysslvrcaches = (nothing, nothing) + odeslvrcache = (reuse, has_explicit, ui_pre, ui, slopes, J, r, sysslvrcaches) + + (odeslvrcache, odeopcache) +end + +function ode_march!( + stateF::NTuple{1,AbstractVector}, + odeslvr::MutableDIMRungeKutta, odeop::ODEOperator, + t0::Real, state0::NTuple{1,AbstractVector}, + odecache +) + # Unpack inputs + u0 = state0[1] + odeslvrcache, odeopcache = odecache + reuse, has_explicit, ui_pre, ui, slopes, J, r, sysslvrcaches = odeslvrcache + sysslvrcache_nl, sysslvrcache_l = sysslvrcaches + + # Unpack solver + sysslvr_nl, sysslvr_l = odeslvr.sysslvr_nl, odeslvr.sysslvr_l + dt, tableau = odeslvr.dt, odeslvr.tableau + A, b, c = get_matrix(tableau), get_weights(tableau), get_nodes(tableau) + + for i in eachindex(c) + # Define scheme + x = slopes[i] + tx = t0 + c[i] * dt + copy!(ui_pre, u0) + for j in 1:i-1 + axpy!(A[i, j] * dt, slopes[j], ui_pre) + end + + # Update ODE operator cache + update_odeopcache!(odeopcache, odeop, tx) + + # Decide whether the stage is explicit or implicit + # The stage becomes explicit when aii = 0 and the operator is quasilinear, + # which is precomputed in has_explicit + aii = A[i, i] + explicit_stage = iszero(aii) && has_explicit + + if explicit_stage + # Define scheme + # Set x to zero to split jacobian and residual + fill!(x, zero(eltype(x))) + usx = (ui_pre, x) + ws = (0, 1) + + # Create and solve stage operator + stageop = LinearStageOperator( + odeop, odeopcache, + tx, usx, ws, + J, r, reuse, sysslvrcache_l + ) + + sysslvrcache_l = solve!(x, sysslvr_l, stageop, sysslvrcache_l) + else + # Define scheme + function usx(x) + copy!(ui, ui_pre) + axpy!(aii * dt, x, ui) + (ui, x) + end + ws = (aii * dt, 1) + + # Create and solve stage operator + stageop = NonlinearStageOperator( + odeop, odeopcache, + tx, usx, ws + ) + + sysslvrcache_nl = solve!(x, sysslvr_nl, stageop, sysslvrcache_nl) + end + end + + # Update state + tF = t0 + dt + stateF = _update_dimrk!(stateF, state0, dt, slopes, b) + + # Pack outputs + sysslvrcaches = (sysslvrcache_nl, sysslvrcache_l) + odeslvrcache = (reuse, has_explicit, ui_pre, ui, slopes, J, r, sysslvrcaches) + odecache = (odeslvrcache, odeopcache) + (tF, stateF, odecache) +end + +############### +# Linear case # +############### +function allocate_odecache( + odeslvr::MutableDIMRungeKutta, odeop::ODEOperator{<:AbstractLinearODE}, + t0::Real, us0::NTuple{1,AbstractVector} +) + u0 = us0[1] + us0N = (u0, u0) + odeopcache = allocate_odeopcache(odeop, t0, us0N) + + ui_pre = zero(u0) + num_stages = length(get_nodes(odeslvr.tableau)) + slopes = [zero(u0) for _ in 1:num_stages] + + stiffness_constant = is_form_constant(odeop, 0) + mass_constant = is_form_constant(odeop, 1) + reuse = (stiffness_constant && mass_constant) + + J = allocate_jacobian(odeop, t0, us0N, odeopcache) + r = allocate_residual(odeop, t0, us0N, odeopcache) + + # Numerical setups for the linear solver + # * If the mass and stiffness matrices are constant, we can reuse numerical + # setups and we allocate one for each distinct aii. + # * Otherwise, there will be no reuse so we only need one numerical setup + # that will be updated. + # To be general, we build a map sysslvrcaches: step -> NumericalSetup. + # We will probably never need more than 256 stages so we can use Int8. + if !reuse + n = 1 + ptrs = fill(Int8(1), num_stages) + else + A = get_matrix(odeslvr.tableau) + d = Dict{eltype(A),Int8}() + n = 0 + ptrs = zeros(Int8, num_stages) + for i in 1:num_stages + aii = A[i, i] + if !haskey(d, aii) + n += 1 + d[aii] = n + end + ptrs[i] = d[aii] + end + end + values = Vector{NumericalSetup}(undef, n) + + sysslvrcaches = CompressedArray(values, ptrs) + odeslvrcache = (reuse, ui_pre, slopes, J, r, sysslvrcaches) + + (odeslvrcache, odeopcache) +end + +function ode_march!( + stateF::NTuple{1,AbstractVector}, + odeslvr::MutableDIMRungeKutta, odeop::ODEOperator{<:AbstractLinearODE}, + t0::Real, state0::NTuple{1,AbstractVector}, + odecache +) + # Unpack inputs + u0 = state0[1] + odeslvrcache, odeopcache = odecache + reuse, ui_pre, slopes, J, r, sysslvrcaches = odeslvrcache + + # Unpack solver + sysslvr = odeslvr.sysslvr_l + dt, tableau = odeslvr.dt, odeslvr.tableau + A, b, c = get_matrix(tableau), get_weights(tableau), get_nodes(tableau) + + for i in eachindex(c) + # Define scheme + # Set x to zero to split jacobian and residual + x = slopes[i] + fill!(x, zero(eltype(x))) + tx = t0 + c[i] * dt + copy!(ui_pre, u0) + for j in 1:i-1 + axpy!(A[i, j] * dt, slopes[j], ui_pre) + end + usx = (ui_pre, x) + ws = (A[i, i] * dt, 1) + + # Update ODE operator cache + update_odeopcache!(odeopcache, odeop, tx) + + # Create and solve stage operator + # sysslvrcaches[i] will be unassigned at the first iteration + sysslvrcache = isassigned(sysslvrcaches, i) ? sysslvrcaches[i] : nothing + stageop = LinearStageOperator( + odeop, odeopcache, + tx, usx, ws, + J, r, reuse, sysslvrcache + ) + + sysslvrcache = solve!(x, sysslvr, stageop, sysslvrcache) + sysslvrcaches = _setindex_all!(sysslvrcaches, sysslvrcache, i) + end + + # Update state + tF = t0 + dt + stateF = _update_dimrk!(stateF, state0, dt, slopes, b) + + # Pack outputs + odeslvrcache = (reuse, ui_pre, slopes, J, r, sysslvrcaches) + odecache = (odeslvrcache, odeopcache) + (tF, stateF, odecache) +end + +######### +# Utils # +######### +function _update_dimrk!( + stateF::NTuple{1,AbstractVector}, state0::NTuple{1,AbstractVector}, + dt::Real, slopes::AbstractVector, b::AbstractVector +) + # uF = u0 + ∑_{1 ≤ i ≤ s} b[i] * dt * slopes[i] + u0 = state0[1] + uF = stateF[1] + + copy!(uF, u0) + for (bi, slopei) in zip(b, slopes) + axpy!(bi * dt, slopei, uF) + end + + (uF,) +end + +function MutableRungeKutta( + sysslvr_nl::NonlinearSolver, sysslvr_l::NonlinearSolver, + dt::Real, tableau::AbstractTableau +) + type = TableauType(tableau) + if type == ExplicitTableau + @notimplemented + elseif type == DiagonallyImplicitTableau + MutableDIMRungeKutta(sysslvr_nl, sysslvr_l, dt, tableau) + elseif type == ImplicitExplicitTableau + @notimplemented + # elseif type == FullyImplicitTableau + # FIMRungeKutta(sysslvr_nl, sysslvr_l, dt, tableau) + end +end + +function MutableRungeKutta( + sysslvr_nl::NonlinearSolver, sysslvr_l::NonlinearSolver, + dt::Real, name::Symbol +) +MutableRungeKutta(sysslvr_nl, sysslvr_l, dt, ButcherTableau(name)) +end + +function MutableRungeKutta(sysslvr_nl::NonlinearSolver, dt::Real, tableau) + MutableRungeKutta(sysslvr_nl, sysslvr_nl, dt, name) +end diff --git a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl new file mode 100644 index 00000000..2c89de8e --- /dev/null +++ b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl @@ -0,0 +1,121 @@ +abstract type StabilisationMethod end + +""" + mutable struct StabilisedReinit{V,M} <: Reinitialiser + +Stabilised FE method for level-set reinitialisation. Artificial viscosity +approach (`ArtificialViscosity`) based on Mallon et al. (2023). DOI: `10.48550/arXiv.2303.13672`. +Interior jump penalty approach (`InteriorPenalty`) adapted from that work and replaces +the artifical viscosity term with an interior jump penalty term. + +# Parameters +- + +""" +mutable struct StabilisedReinit{A,B,C,D,E} <: Reinitialiser + nls::NonlinearSolver + stabilisation_method::A + model::B + space::C + dΩ::D + assembler::Assembler + params::E + cache + function StabilisedReinit(model::B,V_φ::C,dΩ::D,h::Real; + γd = 20.0, + nls = NewtonSolver(LUSolver();maxiter=50,rtol=1.e-14,verbose=true), + assembler=SparseMatrixAssembler(V_φ,V_φ), + stabilisation_method::A = InteriorPenalty(V_φ)) where {A,B,C,D} + params = (;γd,h,order=get_order(V_φ)) + new{A,B,C,D,typeof(params)}(nls,stabilisation_method,model, + V_φ,dΩ,assembler,params,nothing) + end +end + +get_nls(s::StabilisedReinit) = s.nls +get_stabilisation_method(s::StabilisedReinit) = s.stabilisation_method +get_assembler(s::StabilisedReinit) = s.assembler +get_space(s::StabilisedReinit) = s.space +get_model(s::StabilisedReinit) = s.model +get_measure(s::StabilisedReinit) = s.dΩ +get_params(s::StabilisedReinit) = s.params +get_cache(s::StabilisedReinit) = s.cache + +function solve!(s::StabilisedReinit,φh,nls_cache::Nothing) + nls, V_φ, assembler, = s.nls, s.space, s.assembler + # Weak form + res, jac = get_residual_and_jacobian(s,φh) + # Operator and cache + op = get_algebraic_operator(FEOperator(res,jac,V_φ,V_φ,assembler)) + nls_cache = instantiate_caches(get_free_dof_values(φh),nls,op) + s.cache = nls_cache + # Solve + solve!(get_free_dof_values(φh),nls,op,nls_cache) + return φh +end + +function solve!(s::StabilisedReinit,φh,nls_cache) + nls, V_φ, assembler, = s.nls, s.space, s.assembler + # Weak form + res, jac = get_residual_and_jacobian(s,φh) + # Operator and cache + op = get_algebraic_operator(FEOperator(res,jac,V_φ,V_φ,assembler)) + # Solve + solve!(get_free_dof_values(φh),nls,op,nls_cache) + return φh +end + +struct ArtificialViscosity <: StabilisationMethod + stabilisation_coefficent::Number +end + +function get_residual_and_jacobian(s::StabilisedReinit{ArtificialViscosity},φh) + model, dΩ, = s.model, s.dΩ + γd, h, order = s.params + ca = s.stabilisation_method.stabilisation_coefficent + ϵ = 1e-20 + + geo = DiscreteGeometry(φh,model) + cutgeo = cut(model,geo) + Γ = EmbeddedBoundary(cutgeo) + dΓ = Measure(Γ,2*order) + + W(u,∇u) = sign(u) * ∇u / (ϵ + norm(∇u)) + V(w) = ca*h*(sqrt ∘ ( w ⋅ w )) + a(w,u,v) = ∫(v*(W ∘ (w,∇(w))) ⋅ ∇(u) + V(W ∘ (w,∇(w)))*∇(u) ⋅ ∇(v))dΩ + ∫((γd/h)*v*u)dΓ + b(w,v) = ∫((sign ∘ w)*v)dΩ + res(u,v) = a(u,u,v) - b(u,v) + jac(u,du,v) = a(u,du,v) + + return res,jac +end + +struct InteriorPenalty <: StabilisationMethod + dΛ::Measure +end +function InteriorPenalty(V_φ::FESpace) + Λ = SkeletonTriangulation(get_triangulation(V_φ)) + dΛ = Measure(Λ,2get_order(V_φ)) + return InteriorPenalty(dΛ) +end + +function get_residual_and_jacobian(s::StabilisedReinit{InteriorPenalty},φh) + model, dΩ, = s.model, s.dΩ + γd, h, order = s.params + ϵ = 1e-20 + dΛ = s.stabilisation_method.dΛ + + geo = DiscreteGeometry(φh,model) + cutgeo = cut(model,geo) + Γ = EmbeddedBoundary(cutgeo) + dΓ = Measure(Γ,2*order) + + W(u,∇u) = sign(u) * ∇u / (ϵ + norm(∇u)) + V(w) = ca*h*(sqrt ∘ ( w ⋅ w )) + a(w,u,v) = ∫(v*(W ∘ (w,∇(w))) ⋅ ∇(u))dΩ + ∫(h^2*jump(∇(u)) ⋅ jump(∇(v)))dΛ + ∫((γd/h)*v*u)dΓ + b(w,v) = ∫((sign ∘ w)*v)dΩ + res(u,v) = a(u,u,v) - b(u,v) + jac(u,du,v) = a(u,du,v) + + return res,jac +end \ No newline at end of file diff --git a/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl b/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl new file mode 100644 index 00000000..47d279a8 --- /dev/null +++ b/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl @@ -0,0 +1,85 @@ +""" + abstract type Evolver + +Your own unfitted level-set evolution method can be created by implementing +concrete functionality for `solve!`. +""" +abstract type Evolver end + +""" + solve!(::Evolver,φ,args...) + +Evolve the level set function φ according to an Evolution method. +""" +function solve!(::Evolver,φ,args...) + @abstractmethod +end + +""" + abstract type Reinitialiser + +Your own unfitted level-set reinitialisation method can be created by implementing +concrete functionality for `solve!`. +""" +abstract type Reinitialiser end + +""" + solve!(::Reinitialiser,φ,args...) + +Reinitialise the level set function φ according to an Reinitialiser method. +""" +function solve!(::Reinitialiser,φ,args...) + @abstractmethod +end + +""" + struct UnfittedFEEvolution{A<:Evolver,B<:Reinitialiser} <: LevelSetEvolution + +Wrapper for unfitted evolution and reinitialisation methods. +""" +struct UnfittedFEEvolution{A<:Evolver,B<:Reinitialiser} <: LevelSetEvolution + evolver::A + reinitialiser::B +end + +function evolve!(s::UnfittedFEEvolution,φ::T,vel::M,γ,args...) where + {T<:AbstractVector,M<:AbstractVector} + V_φ = get_space(s.evolver) + φh = FEFunction(V_φ,φ) + velh = FEFunction(V_φ,vel) + evolve!(s,φh,velh,γ,args...) +end + +function evolve!(s::UnfittedFEEvolution,φ::T,velh,γ,args...) where T<:AbstractVector +V_φ = get_space(s.evolver) +φh = FEFunction(V_φ,φ) +evolve!(s,φh,velh,γ,args...) +end + +function evolve!(s::UnfittedFEEvolution,φh,velh,γ,args...) + cache = get_cache(s.evolver) + solve!(s.evolver,φh,velh,γ,cache) + return get_free_dof_values(φh) +end + +function reinit!(s::UnfittedFEEvolution,φ::T,args...) where T<:AbstractVector + V_φ = get_space(s.reinitialiser) + φh = FEFunction(V_φ,φ) + reinit!(s,φh,args...) +end + +function reinit!(s::UnfittedFEEvolution,φh,args...) + cache = get_cache(s.reinitialiser) + solve!(s.reinitialiser,φh,cache) + return get_free_dof_values(φh) +end + +function get_dof_Δ(s::UnfittedFEEvolution) + V_φ = get_space(s.evolver) + _,h,_ = get_params(s.evolver) + return h/get_order(V_φ) +end + +include("CutFEMEvolve.jl") +include("StabilisedReinit.jl") +include("MutableRungeKutta.jl") \ No newline at end of file diff --git a/src/Utilities.jl b/src/Utilities.jl index df37fa35..a3bc5786 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -231,4 +231,14 @@ function _zerocrossing_period(v;t=0:length(v)-1,line=sum(v)/length(v)) inds = findall(@. ≥(line, $@view(v[2:end])) & <(line, $@view(v[1:end-1]))) difft = diff(t[inds]) sum(difft)/length(difft) +end + +import Gridap.ReferenceFEs: get_order +function get_order(space::FESpace) + return get_order(first(Gridap.CellData.get_data(get_fe_basis(space)))) +end + +function get_order(space::DistributedFESpace) + order = map(get_order,local_views(space)) + return getany(order) end \ No newline at end of file