Added UnfittedEvolution to src and tests

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

@@ -1,108 +1,108 @@
-# 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)]
-)
-
-
-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.
-
+# 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)]
+)
+
+
+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.
+

src/GridapTopOpt.jl

@@ -16,6 +16,7 @@ using Gridap.ReferenceFEs, Gridap.FESpaces,  Gridap.MultiField
 
 using Gridap.Geometry: get_faces, num_nodes
 using Gridap.FESpaces: get_assembly_strategy
+using Gridap.ODEs: ODESolver
 using Gridap: writevtk
 
 using GridapDistributed
@@ -38,6 +39,7 @@ using GridapEmbedded.Interfaces: SubFacetData, SubCellTriangulation, SubFacetTri
 using JLD2: save_object, load_object, jldsave
 
 import Base: +
+import Gridap: solve!
 
 include("GridapExtensions.jl")
 
@@ -65,6 +67,11 @@ export isotropic_elast_tensor
 include("LevelSetEvolution/LevelSetEvolution.jl")
 export HamiltonJacobiEvolution
 export FirstOrderStencil
+export UnfittedFEEvolution
+export CutFEMEvolve
+export StabilisedReinit
+export ArtificialViscosity
+export InteriorPenalty
 export evolve!
 export reinit!