diff --git a/scripts/Embedded/Bugs/connor_issue.jl b/scripts/Embedded/Bugs/connor_issue.jl new file mode 100644 index 0000000..566acaa --- /dev/null +++ b/scripts/Embedded/Bugs/connor_issue.jl @@ -0,0 +1,54 @@ +module PoissonMultiFieldTests + +using Gridap +using Test + +u1(x) = x[1]^2 +u2(x) = x[2]*x[1] + +f1(x) = -Δ(u1)(x) +f2(x) = u1(x)*-(Δ(u2)(x)) + +domain = (0,1,0,1) +cells = (2,2) +model = CartesianDiscreteModel(domain,cells) + +order = 2 +V = FESpace(model, ReferenceFE(lagrangian,Float64,order),conformity=:H1,dirichlet_tags="boundary") +U1 = TrialFESpace(V,u1) +U2 = TrialFESpace(V,u2) + +Ω = Triangulation(model) +Γ = BoundaryTriangulation(model,tags=1) +n_Γ = get_normal_vector(Γ) + +degree = 2*order +dΩ = Measure(Ω,degree) +dΓ = Measure(Γ,degree) + +a1(u1,v1) = ∫( ∇(v1)⋅∇(u1) )*dΩ +l1(v1) = ∫( v1*f1 )*dΩ + +op1 = AffineFEOperator(a1,l1,U1,V) +u1h = solve(op1) + +a2(u2,v2) = ∫( u1h * ∇(v2)⋅∇(u2) )*dΩ +l2(v2) = ∫( v2*f2 )*dΩ + +op2 = AffineFEOperator(a2,l2,U2,V) +u2h = solve(op2) + +a3((u1,u2),(v1,v2)) = ∫( ∇(v1)⋅∇(u1) )dΩ + ∫( u1*∇(v2)⋅∇(u2))dΩ +l3((v1,v2)) = ∫( v1*f1 )*dΩ + ∫( v2*f2 )*dΩ + +X = MultiFieldFESpace([U1,U2]) +Y = MultiFieldFESpace([V,V]) + +op3 = FEOperator((u,v)->a3(u,v)-l3(v),X,Y) +solver = NLSolver(LUSolver(), show_trace=true, method=:newton, iterations=5) +(u1hm,u2hm) = solve(solver,op3) + +op3 = AffineFEOperator(a3,l3,X,Y) +solve(op3) + +end \ No newline at end of file diff --git a/scripts/Embedded/Examples/CutFEM_2d_thermal_(island_checking).jl b/scripts/Embedded/Examples/CutFEM_2d_thermal_(island_checking).jl index f6f9ddd..f7dc667 100644 --- a/scripts/Embedded/Examples/CutFEM_2d_thermal_(island_checking).jl +++ b/scripts/Embedded/Examples/CutFEM_2d_thermal_(island_checking).jl @@ -68,7 +68,7 @@ a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + l(v,φ) = ∫(v)dΓ_N ## Optimisation functionals -J(u,φ) = a(u,u,φ) +J(u,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ dVol(q,u,φ) = ∫(-1/vol_D*q/(1e-20 + norm ∘ (∇(φ))))Ωs.dΓ diff --git a/scripts/Embedded/Examples/FCM_2d_thermal.jl b/scripts/Embedded/Examples/FCM_2d_thermal.jl index 8dfd54a..eff48be 100644 --- a/scripts/Embedded/Examples/FCM_2d_thermal.jl +++ b/scripts/Embedded/Examples/FCM_2d_thermal.jl @@ -63,12 +63,12 @@ V_φ = TestFESpace(model,reffe_scalar) end ## Weak form -const ϵ = 1e-3 +ϵ = 1e-3 a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + ∫(ϵ*∇(v)⋅∇(u))Ωs.dΩout l(v,φ) = ∫(v)dΓ_N ## Optimisation functionals -J(u,φ) = a(u,u,φ) +J(u,φ) = ∫(∇(u)⋅∇(u))Ωs.dΩin Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ dVol(q,u,φ) = ∫(-1/vol_D*q/(norm ∘ (∇(φ))))Ωs.dΓ diff --git a/scripts/Embedded/Examples/FCM_2d_thermal_MPI.jl b/scripts/Embedded/Examples/FCM_2d_thermal_MPI.jl index c1a9d37..86f2b10 100644 --- a/scripts/Embedded/Examples/FCM_2d_thermal_MPI.jl +++ b/scripts/Embedded/Examples/FCM_2d_thermal_MPI.jl @@ -77,7 +77,7 @@ function main(mesh_partition,distribute,el_size,path) l(v,φ) = ∫(v)dΓ_N ## Optimisation functionals - J(u,φ) = a(u,u,φ) + J(u,φ) = ∫(∇(u)⋅∇(u))Ωs.dΩin Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ dVol(q,u,φ) = ∫(-1/vol_D*q/(norm ∘ (∇(φ))))Ωs.dΓ diff --git a/scripts/Embedded/Examples/FCM_2d_thermal_with_NonDesignDom.jl b/scripts/Embedded/Examples/FCM_2d_thermal_with_NonDesignDom.jl index d4c7162..f75154d 100644 --- a/scripts/Embedded/Examples/FCM_2d_thermal_with_NonDesignDom.jl +++ b/scripts/Embedded/Examples/FCM_2d_thermal_with_NonDesignDom.jl @@ -70,7 +70,7 @@ a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + ∫(ϵ*∇(v)⋅∇(u))Ωs.dΩout l(v,φ) = ∫(v)dΓ_N ## Optimisation functionals -J(u,φ) = a(u,u,φ) +J(u,φ) = ∫(∇(u)⋅∇(u))Ωs.dΩin Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ dVol(q,u,φ) = ∫(-1/vol_D*q/(norm ∘ (∇(φ))))Ωs.dΓ diff --git a/scripts/Embedded/Examples/fsi/TO-2-CutFEM-stokes_fsi.jl b/scripts/Embedded/Examples/fsi/TO-2-CutFEM-stokes_fsi.jl new file mode 100644 index 0000000..f26f33e --- /dev/null +++ b/scripts/Embedded/Examples/fsi/TO-2-CutFEM-stokes_fsi.jl @@ -0,0 +1,232 @@ +using Gridap, Gridap.Geometry, Gridap.Adaptivity +using GridapEmbedded, GridapEmbedded.LevelSetCutters +using GridapTopOpt + +using GridapTopOpt: StateParamIntegrandWithMeasure + +path = "./results/TO-2-CutFEM-stokes_fsi/" +mkpath(path) + +n = 100 +γ_evo = 0.05 +max_steps = floor(Int,n/5) +vf = 0.03 +α_coeff = 2 +iter_mod = 1 + +# Cut the background model +L = 2.0 +H = 0.5 +x0 = 0.5 +l = 0.4 +w = 0.05 +a = 0.3 +b = 0.03 +vol_D = L*H + +L - l/2 + +partition = (4n,n) +D = length(partition) +_model = CartesianDiscreteModel((0,L,0,H),partition) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model + +el_Δ = get_el_Δ(_model) +h = maximum(el_Δ) +h_refine = maximum(el_Δ)/2 + +f_Γ_D(x) = x[1] ≈ 0 +f_Γ_N(x) = x[1] ≈ L +f_Γ_NoSlipTop(x) = x[2] ≈ H +f_Γ_NoSlipBottom(x) = x[2] ≈ 0 +f_NonDesign(x) = ((x0 - w/2 - eps() <= x[1] <= x0 + w/2 + eps() && 0.0 <= x[2] <= a + eps()) || + (x0 - l/2 - eps() <= x[1] <= x0 + l/2 + eps() && 0.0 <= x[2] <= b + eps())) +update_labels!(1,model,f_Γ_D,"Gamma_f_D") +update_labels!(2,model,f_Γ_N,"Gamma_f_N") +update_labels!(3,model,f_Γ_NoSlipTop,"Gamma_NoSlipTop") +update_labels!(4,model,f_Γ_NoSlipBottom,"Gamma_NoSlipBottom") +update_labels!(5,model,f_NonDesign,"NonDesign") +update_labels!(6,model,x->f_NonDesign(x) && f_Γ_NoSlipBottom(x),"Gamma_s_D") +# writevtk(model,path*"model") + +# Cut the background model +reffe_scalar = ReferenceFE(lagrangian,Float64,1) +V_φ = TestFESpace(model,reffe_scalar) +V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["NonDesign","Gamma_s_D"]) +U_reg = TrialFESpace(V_reg) + +f0((x,y),W,H) = max(2/W*abs(x-x0),1/(H/2+1)*abs(y-H/2+1))-1 +f1((x,y),q,r) = - cos(q*π*x)*cos(q*π*y)/q - r/q +fin(x) = f0(x,l,a) +fsolid(x) = min(f0(x,l,b),f0(x,w,a)) +fholes((x,y),q,r) = max(f1((x,y),q,r),f1((x-1/q,y),q,r)) +φf(x) = min(max(fin(x),fholes(x,15,0.5)),fsolid(x)) +φh = interpolate(φf,V_φ) +# writevtk(get_triangulation(φh),path*"initial_lsf",cellfields=["φ"=>φh]) + +# Setup integration meshes and measures +order = 1 +degree = 2*order + +Ω_bg = Triangulation(model) +dΩ_bg = Measure(Ω_bg,2*order) +Γf_D = BoundaryTriangulation(model,tags="Gamma_f_D") +Γf_N = BoundaryTriangulation(model,tags="Gamma_f_N") +dΓf_D = Measure(Γf_D,degree) +dΓf_N = Measure(Γf_N,degree) +Ω = EmbeddedCollection(model,φh) do cutgeo,_ + Ωs = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL),V_φ) + Ωf = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ) + Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ) + Γg = GhostSkeleton(cutgeo) + (; + :Ωs => Ωs, + :dΩs => Measure(Ωs,degree), + :Ωf => Ωf, + :dΩf => Measure(Ωf,degree), + :Γg => Γg, + :dΓg => Measure(Γg,degree), + :n_Γg => get_normal_vector(Γg), + :Γ => Γ, + :dΓ => Measure(Γ,degree), + :n_Γ => get_normal_vector(Γ.trian), + :Ω_act_s => Triangulation(cutgeo,ACTIVE), + :Ω_act_f => Triangulation(cutgeo,ACTIVE_OUT), + :χ_s => GridapTopOpt.get_isolated_volumes_mask(cutgeo,["Gamma_s_D"];IN_is=IN), + :χ_f => GridapTopOpt.get_isolated_volumes_mask(cutgeo,["Gamma_f_D"];IN_is=OUT) + ) +end + +# Setup spaces +uin(x) = VectorValue(16x[2]*(H-x[2]),0.0) + +reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P) +reffe_p = ReferenceFE(lagrangian,Float64,order,space=:P) +reffe_d = ReferenceFE(lagrangian,VectorValue{D,Float64},order) + +function build_spaces(Ω_act_s,Ω_act_f) + # Test spaces + V = TestFESpace(Ω_act_f,reffe_u,conformity=:H1, + dirichlet_tags=["Gamma_f_D","Gamma_NoSlipTop","Gamma_NoSlipBottom","Gamma_s_D"]) + Q = TestFESpace(Ω_act_f,reffe_p,conformity=:H1) + T = TestFESpace(Ω_act_s,reffe_d,conformity=:H1,dirichlet_tags=["Gamma_s_D"]) + + # Trial spaces + U = TrialFESpace(V,[uin,VectorValue(0.0,0.0),VectorValue(0.0,0.0),VectorValue(0.0,0.0)]) + P = TrialFESpace(Q) + R = TrialFESpace(T) + + # Multifield spaces + X = MultiFieldFESpace([U,P,R]) + Y = MultiFieldFESpace([V,Q,T]) + return X,Y +end + +# Weak form +## Fluid +# Properties +Re = 60 # Reynolds number +ρ = 1.0 # Density +cl = H # Characteristic length +u0_max = 1.0 # Maximum velocity on inlet +μ = 1.0#ρ*cl*u0_max/Re # Viscosity +# Stabilization parameters +β0 = 0.25 +β1 = 0.2 +β2 = 0.1 +β3 = 0.05 +γ = 100.0 + +# Terms +σf_n(u,p,n) = μ*∇(u)⋅n - p*n +a_Ω(u,v) = μ*(∇(u) ⊙ ∇(v)) +b_Ω(v,p) = - (∇⋅v)*p +c_Γi(p,q) = (β0*h)*jump(p)*jump(q) # this will vanish for p∈P1 +c_Ω(p,q) = (β1*h^2)*∇(p)⋅∇(q) +a_Γ(u,v) = - (Ω.n_Γ⋅∇(u))⋅v - u⋅(Ω.n_Γ⋅∇(v)) + (γ/h)*u⋅v +b_Γ(v,p) = (Ω.n_Γ⋅v)*p +i_Γg(u,v) = (β2*h)*jump(Ω.n_Γg⋅∇(u))⋅jump(Ω.n_Γg⋅∇(v)) +j_Γg(p,q) = (β3*h^3)*jump(Ω.n_Γg⋅∇(p))*jump(Ω.n_Γg⋅∇(q)) + c_Γi(p,q) + +a_fluid((u,p),(v,q)) = + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p)-c_Ω(p,q) )Ω.dΩf + + ∫( a_Γ(u,v)+b_Γ(u,q)+b_Γ(v,p) )Ω.dΓ + + ∫( i_Γg(u,v) - j_Γg(p,q) )Ω.dΓg + + ∫(Ω.χ_f*(p*q + u⋅v))Ω.dΩf # Isolated volume term + +## Structure +# Stabilization and material parameters +function lame_parameters(E,ν) + λ = (E*ν)/((1+ν)*(1-2*ν)) + μ = E/(2*(1+ν)) + (λ, μ) +end +λs, μs = lame_parameters(1.0,0.3) #0.1,0.05) +γg = (λs + 2μs)*0.1 +# Terms +σ(ε) = λs*tr(ε)*one(ε) + 2*μs*ε +a_solid(d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))Ω.dΩs + + ∫((γg*h^3)*jump(Ω.n_Γg⋅∇(s)) ⋅ jump(Ω.n_Γg⋅∇(d)))Ω.dΓg + + ∫(Ω.χ_s*d⋅s)Ω.dΩs # Isolated volume term + +## Full problem +# minus sign because of the normal direction +function a_coupled((u,p,d),(v,q,s),φ) + n_AD = get_normal_vector(Ω.Γ) + return a_fluid((u,p),(v,q)) + a_solid(d,s) + ∫(-σf_n(u,p,n_AD) ⋅ s)Ω.dΓ +end +l_coupled((v,q,s),φ) = ∫(0.0q)Ω.dΓ + +## Optimisation functionals +J_pres((u,p,d),φ) = ∫(p)dΓf_D - ∫(p)dΓf_N +J_comp((u,p,d),φ) = ∫(ε(d) ⊙ (σ ∘ ε(d)))Ω.dΩs +Vol((u,p,d),φ) = ∫(100*1/vol_D)Ω.dΩs - ∫(100*vf/vol_D)dΩ_bg + +## Setup solver and FE operators +state_collection = EmbeddedCollection(model,φh) do _,_ + X,Y = build_spaces(Ω.Ω_act_s,Ω.Ω_act_f) + state_map = AffineFEStateMap(a_coupled,l_coupled,X,Y,V_φ,U_reg,φh) + (; + :state_map => state_map, + :J => StateParamIntegrandWithMeasure(J_comp,state_map), + :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,]) + ) +end +pcfs = EmbeddedPDEConstrainedFunctionals(state_collection) + +## Evolution Method +evo = CutFEMEvolve(V_φ,Ω,dΩ_act,h;max_steps) +reinit = StabilisedReinit(V_φ,Ω,dΩ_act,h;stabilisation_method=ArtificialViscosity(3.0)) +ls_evo = UnfittedFEEvolution(evo,reinit) + +## Hilbertian extension-regularisation problems +α = α_coeff*h_refine +a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ_act; +vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) + +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; + γ=γ_evo,verbose=true,constraint_names=[:Vol]) +for (it,(uh,ph,dh),φh) in optimiser + if iszero(it % iter_mod) + writevtk(Ω_bg,path*"Omega_act_$it", + cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"ph"=>ph, + "dh"=>dh,"χ_s"=>Ω.χ_s,"χ_f"=>Ω.χ_f]) + writevtk(Ω.Ωf,path*"Omega_f_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) + writevtk(Ω.Ωs,path*"Omega_s_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) + writevtk(Ω.Γ,path*"Gamma_$it",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅Ω.n_Γ,"σf_n"=>σf_n(uh,ph,φh)]) + error() + end + write_history(path*"/history.txt",optimiser.history) +end +it = get_history(optimiser).niter; uh,ph,dh = get_state(pcfs) +writevtk(Ω_bg,path*"Omega_act_$it", + cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ω.Ωf,path*"Omega_f_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ω.Ωs,path*"Omega_s_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ω.Γ,path*"Gamma_$it",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅Ω.n_Γ,"σf_n"=>σf_n(uh,ph,φh)]) \ No newline at end of file diff --git a/scripts/Embedded/Examples/fsi/TO-5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl b/scripts/Embedded/Examples/fsi/TO-5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl new file mode 100644 index 0000000..99d9707 --- /dev/null +++ b/scripts/Embedded/Examples/fsi/TO-5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl @@ -0,0 +1,204 @@ +using Gridap, Gridap.Geometry, Gridap.Adaptivity +using GridapEmbedded, GridapEmbedded.LevelSetCutters +using GridapTopOpt + +path = "./results/TO-6-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi/" +mkpath(path) + +n = 100 +γ_evo = 0.05 +max_steps = floor(Int,n/5) +vf = 0.03 +α_coeff = 2 +iter_mod = 1 + +# Cut the background model +L = 2.0 +H = 0.5 +x0 = 0.5 +l = 0.4 +w = 0.05 +a = 0.3 +b = 0.03 + +partition = (4n,n) +D = length(partition) +_model = CartesianDiscreteModel((0,L,0,H),partition) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model + +el_Δ = get_el_Δ(_model) +h = maximum(el_Δ) +h_refine = maximum(el_Δ)/2 + +f_Γ_D(x) = x[1] ≈ 0 +f_Γ_N(x) = x[1] ≈ L +f_Γ_NoSlipTop(x) = x[2] ≈ H +f_Γ_NoSlipBottom(x) = x[2] ≈ 0 +f_NonDesign(x) = ((x0 - w/2 - eps() <= x[1] <= x0 + w/2 + eps() && 0.0 <= x[2] <= a + eps()) || + (x0 - l/2 - eps() <= x[1] <= x0 + l/2 + eps() && 0.0 <= x[2] <= b + eps())) +update_labels!(1,model,f_Γ_D,"Gamma_f_D") +update_labels!(2,model,f_Γ_N,"Gamma_f_N") +update_labels!(3,model,f_Γ_NoSlipTop,"Gamma_NoSlipTop") +update_labels!(4,model,f_Γ_NoSlipBottom,"Gamma_NoSlipBottom") +update_labels!(5,model,f_NonDesign,"NonDesign") +update_labels!(6,model,x->f_NonDesign(x) && f_Γ_NoSlipBottom(x),"Gamma_s_D") +# writevtk(model,path*"model") + +# Cut the background model +reffe_scalar = ReferenceFE(lagrangian,Float64,1) +V_φ = TestFESpace(model,reffe_scalar) +V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["NonDesign","Gamma_s_D"]) +U_reg = TrialFESpace(V_reg) + +f0((x,y),W,H) = max(2/W*abs(x-x0),1/(H/2+1)*abs(y-H/2+1))-1 +f1((x,y),q,r) = - cos(q*π*x)*cos(q*π*y)/q - r/q +fin(x) = f0(x,l,a) +fsolid(x) = min(f0(x,l,b),f0(x,w,a)) +fholes((x,y),q,r) = max(f1((x,y),q,r),f1((x-1/q,y),q,r)) +φf(x) = min(max(fin(x),fholes(x,15,0.5)),fsolid(x)) +φh = interpolate(φf,V_φ) +# writevtk(get_triangulation(φh),path*"initial_lsf",cellfields=["φ"=>φh]) + +# Setup integration meshes and measures +order = 1 +degree = 2*order + +Ω_act = Triangulation(model) +dΩ_act = Measure(Ω_act,degree) +vol_D = L*H + +Γf_D = BoundaryTriangulation(model,tags="Gamma_f_D") +Γf_N = BoundaryTriangulation(model,tags="Gamma_f_N") +dΓf_D = Measure(Γf_D,degree) +dΓf_N = Measure(Γf_N,degree) +Ω = EmbeddedCollection(model,φh) do cutgeo,_ + Ωs = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL),V_φ) + Ωf = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ) + Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ) + Γg = GhostSkeleton(cutgeo) + Ωact = Triangulation(cutgeo,ACTIVE) + (; + :Ωs => Ωs, + :dΩs => Measure(Ωs,degree), + :Ωf => Ωf, + :dΩf => Measure(Ωf,degree), + :Γg => Γg, + :dΓg => Measure(Γg,degree), + :n_Γg => get_normal_vector(Γg), + :Γ => Γ, + :dΓ => Measure(Γ,degree), + :n_Γ => get_normal_vector(Γ.trian), + :Ωact => Ωact + ) +end + +# Setup FESpace +uin(x) = VectorValue(16x[2]*(H-x[2]),0.0) + +reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P) +reffe_p = ReferenceFE(lagrangian,Float64,order,space=:P) +reffe_d = ReferenceFE(lagrangian,VectorValue{D,Float64},order) + +V = TestFESpace(Ω_act,reffe_u,conformity=:H1, + dirichlet_tags=["Gamma_f_D","Gamma_NoSlipTop","Gamma_NoSlipBottom","Gamma_s_D"]) +Q = TestFESpace(Ω_act,reffe_p,conformity=:H1) +T = TestFESpace(Ω_act ,reffe_d,conformity=:H1,dirichlet_tags=["Gamma_s_D"]) + +U = TrialFESpace(V,[uin,VectorValue(0.0,0.0),VectorValue(0.0,0.0),VectorValue(0.0,0.0)]) +P = TrialFESpace(Q) +R = TrialFESpace(T) + +X = MultiFieldFESpace([U,P,R]) +Y = MultiFieldFESpace([V,Q,T]) + +# Weak form +## Fluid +# Properties +Re = 60 # Reynolds number +ρ = 1.0 # Density +cl = H # Characteristic length +u0_max = maximum(abs,get_dirichlet_dof_values(U)) +μ = 1.0#ρ*cl*u0_max/Re # Viscosity +# Stabilization parameters +# TODO: Need to look at this much closer +γ = 10^16#10.0/μ +β1 = 0.2#1/μ + +# Terms +σf_n(u,p,n) = μ*∇(u) ⋅ n - p*n +a_Ω(u,v) = μ*(∇(u) ⊙ ∇(v)) +b_Ω(v,p) = - (∇⋅v)*p +c_Ω(p,q) = (β1*h^2)*∇(p)⋅∇(q) + +a_fluid((u,p),(v,q)) = + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p)-c_Ω(p,q) ) * Ω.dΩf + + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p)-c_Ω(p,q)+(γ/h)*u⋅v ) * Ω.dΩs + +## Structure +# Stabilization and material parameters +function lame_parameters(E,ν) + λ = (E*ν)/((1+ν)*(1-2*ν)) + μ = E/(2*(1+ν)) + (λ, μ) +end +λs, μs = lame_parameters(1.0,0.3) #0.1,0.05) +ϵ = (λs + 2μs)*1e-3 +# Terms +σ(ε) = λs*tr(ε)*one(ε) + 2*μs*ε +a_solid(d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))Ω.dΩs + + ∫(ϵ*(ε(s) ⊙ (σ ∘ ε(d))))Ω.dΩf # Ersatz + +## Full problem +vec0 = VectorValue(0.0,0.0) +# minus sign because of the normal direction +function a_coupled((u,p,d),(v,q,s),φ) + n_AD = get_normal_vector(Ω.Γ) + return a_fluid((u,p),(v,q)) + a_solid(d,s) + + ∫(-σf_n(u,p,n_AD) ⋅ s)Ω.dΓ + ∫((-σf_n(u,p,vec0) ⋅ s))dΩ_act +end +l_coupled((v,q,s),φ) = ∫(0.0q)dΩ_act + +## Optimisation functionals +J_pres((u,p,d),φ) = ∫(p)dΓf_D - ∫(p)dΓf_N +J_comp((u,p,d),φ) = ∫(ε(d) ⊙ (σ ∘ ε(d)))Ω.dΩs +Vol((u,p,d),φ) = ∫(100*1/vol_D)Ω.dΩs - ∫(100*vf/vol_D)dΩ_act + +## Setup solver and FE operators +state_map = AffineFEStateMap(a_coupled,l_coupled,X,Y,V_φ,U_reg,φh) +pcfs = PDEConstrainedFunctionals(J_pres,[Vol],state_map) + +## Evolution Method +evo = CutFEMEvolve(V_φ,Ω,dΩ_act,h;max_steps) +reinit = StabilisedReinit(V_φ,Ω,dΩ_act,h;stabilisation_method=ArtificialViscosity(3.0)) +ls_evo = UnfittedFEEvolution(evo,reinit) + +## Hilbertian extension-regularisation problems +α = α_coeff*h_refine +a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ_act; +vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) + +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; + γ=γ_evo,verbose=true,constraint_names=[:Vol]) +for (it,(uh,ph,dh),φh) in optimiser + if iszero(it % iter_mod) + writevtk(Ω_act,path*"Omega_act_$it", + cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"ph"=>ph,"dh"=>dh]) + writevtk(Ω.Ωf,path*"Omega_f_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) + writevtk(Ω.Ωs,path*"Omega_s_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) + writevtk(Ω.Γ,path*"Gamma_$it",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅Ω.n_Γ,"σf_n"=>σf_n(uh,ph,φh)]) + error() + end + write_history(path*"/history.txt",optimiser.history) +end +it = get_history(optimiser).niter; uh,ph,dh = get_state(pcfs) +writevtk(Ω_act,path*"Omega_act_$it", + cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ω.Ωf,path*"Omega_f_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ω.Ωs,path*"Omega_s_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ω.Γ,path*"Gamma_$it",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅Ω.n_Γ,"σf_n"=>σf_n(uh,ph,φh)]) \ No newline at end of file diff --git a/scripts/Embedded/Examples/fsi/TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl b/scripts/Embedded/Examples/fsi/TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl index 645c2be..44ede39 100644 --- a/scripts/Embedded/Examples/fsi/TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl +++ b/scripts/Embedded/Examples/fsi/TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl @@ -2,20 +2,26 @@ using Gridap, Gridap.Geometry, Gridap.Adaptivity using GridapEmbedded, GridapEmbedded.LevelSetCutters using GridapTopOpt -path = "./results/fsi testing/" +path = "./results/TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi/" mkpath(path) +n = 100 +γ_evo = 0.1 +max_steps = floor(Int,n/5) +vf = 0.03 +α_coeff = 2 +iter_mod = 1 + # Cut the background model -L = 1.0#2.0 -H = 1.0#0.5 +L = 2.0 +H = 0.5 x0 = 0.5 l = 0.4 -w = 0.1 -a = 0.5#0.3 -b = 0.01 +w = 0.05 +a = 0.3 +b = 0.03 -n = 100 -partition = (n,n) +partition = (4n,n) D = length(partition) _model = CartesianDiscreteModel((0,L,0,H),partition) base_model = UnstructuredDiscreteModel(_model) @@ -24,61 +30,83 @@ model = ref_model.model el_Δ = get_el_Δ(_model) h = maximum(el_Δ) +h_refine = maximum(el_Δ)/2 + f_Γ_D(x) = x[1] ≈ 0 +f_Γ_N(x) = x[1] ≈ L f_Γ_NoSlipTop(x) = x[2] ≈ H f_Γ_NoSlipBottom(x) = x[2] ≈ 0 -f_NonDesign(x) = ((x0 - w/2 <= x[1] <= x0 + w/2 && 0.0 <= x[2] <= a) && - x0 - w/2 <= x[1] <= x0 + w/2 && 0.0 <= x[2] <= b) -update_labels!(1,model,f_Γ_D,"Gamma_D") -update_labels!(2,model,f_Γ_NoSlipTop,"Gamma_NoSlipTop") -update_labels!(3,model,f_Γ_NoSlipBottom,"Gamma_NoSlipBottom") -update_labels!(4,model,f_NonDesign,"NonDesign") +f_NonDesign(x) = ((x0 - w/2 - eps() <= x[1] <= x0 + w/2 + eps() && 0.0 <= x[2] <= a + eps()) || + (x0 - l/2 - eps() <= x[1] <= x0 + l/2 + eps() && 0.0 <= x[2] <= b + eps())) +update_labels!(1,model,f_Γ_D,"Gamma_f_D") +update_labels!(2,model,f_Γ_N,"Gamma_f_N") +update_labels!(3,model,f_Γ_NoSlipTop,"Gamma_NoSlipTop") +update_labels!(4,model,f_Γ_NoSlipBottom,"Gamma_NoSlipBottom") +update_labels!(5,model,f_NonDesign,"NonDesign") +update_labels!(6,model,x->f_NonDesign(x) && f_Γ_NoSlipBottom(x),"Gamma_s_D") +# writevtk(model,path*"model") # Cut the background model reffe_scalar = ReferenceFE(lagrangian,Float64,1) V_φ = TestFESpace(model,reffe_scalar) -φh = interpolate(x->-max(20*abs(x[1]-0.5),3*abs(x[2]-0.2))+1,V_φ) -geo = DiscreteGeometry(φh,model) -cutgeo = cut(model,geo) -cutgeo_facets = cut_facets(model,geo) - -# Generate the "active" model -Ω_act = Triangulation(model) - -# Setup integration meshes -Ω = Triangulation(cutgeo,PHYSICAL) -Ωout = Triangulation(cutgeo,PHYSICAL_OUT) -Γ = EmbeddedBoundary(cutgeo) -Γg = GhostSkeleton(cutgeo) -Γi = SkeletonTriangulation(cutgeo_facets) - -# Setup normal vectors -n_Γ = get_normal_vector(Γ) -n_Γg = get_normal_vector(Γg) -n_Γi = get_normal_vector(Γi) - -# Setup Lebesgue measures +V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["NonDesign","Gamma_s_D"]) +U_reg = TrialFESpace(V_reg) + +f0((x,y),W,H) = max(2/W*abs(x-x0),1/(H/2+1)*abs(y-H/2+1))-1 +f1((x,y),q,r) = - cos(q*π*x)*cos(q*π*y)/q - r/q +fin(x) = f0(x,l,a) +fsolid(x) = min(f0(x,l,b),f0(x,w,a)) +fholes((x,y),q,r) = max(f1((x,y),q,r),f1((x-1/q,y),q,r)) +φf(x) = min(max(fin(x),fholes(x,15,0.5)),fsolid(x)) +φh = interpolate(φf,V_φ) +# writevtk(get_triangulation(φh),path*"initial_lsf",cellfields=["φ"=>φh]) + +# Setup integration meshes and measures order = 2 degree = 2*order -dΩ = Measure(Ω,degree) -dΩout = Measure(Ωout,degree) -dΓ = Measure(Γ,degree) -dΓg = Measure(Γg,degree) -dΓi = Measure(Γi,degree) -# Setup FESpace +Ω_act = Triangulation(model) +dΩ_act = Measure(Ω_act,degree) +vol_D = L*H + +Γf_D = BoundaryTriangulation(model,tags="Gamma_f_D") +Γf_N = BoundaryTriangulation(model,tags="Gamma_f_N") +dΓf_D = Measure(Γf_D,degree) +dΓf_N = Measure(Γf_N,degree) +Ω = EmbeddedCollection(model,φh) do cutgeo,_ + Ωs = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL),V_φ) + Ωf = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ) + Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ) + Γg = GhostSkeleton(cutgeo) + Ωact = Triangulation(cutgeo,ACTIVE) + (; + :Ωs => Ωs, + :dΩs => Measure(Ωs,degree), + :Ωf => Ωf, + :dΩf => Measure(Ωf,degree), + :Γg => Γg, + :dΓg => Measure(Γg,degree), + :n_Γg => get_normal_vector(Γg), + :Γ => Γ, + :dΓ => Measure(Γ,degree), + :n_Γ => get_normal_vector(Γ.trian), + :Ωact => Ωact + ) +end -uin(x) = VectorValue(x[2]*(1-x[2]),0.0) +# Setup FESpace +uin(x) = VectorValue(16x[2]*(H-x[2]),0.0) reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P) reffe_p = ReferenceFE(lagrangian,Float64,order-1,space=:P) -reffe_d = ReferenceFE(lagrangian,VectorValue{D,Float64},order) +reffe_d = ReferenceFE(lagrangian,VectorValue{D,Float64},order-1) -V = TestFESpace(Ω_act,reffe_u,conformity=:H1,dirichlet_tags=["Gamma_D","Gamma_NoSlipTop","Gamma_NoSlipBottom"]) +V = TestFESpace(Ω_act,reffe_u,conformity=:H1, + dirichlet_tags=["Gamma_f_D","Gamma_NoSlipTop","Gamma_NoSlipBottom","Gamma_s_D"]) Q = TestFESpace(Ω_act,reffe_p,conformity=:C0) -T = TestFESpace(Ω_act ,reffe_d,conformity=:H1,dirichlet_tags=["Gamma_NoSlipBottom"]) +T = TestFESpace(Ω_act ,reffe_d,conformity=:H1,dirichlet_tags=["Gamma_s_D"]) -U = TrialFESpace(V,[uin,VectorValue(0.0,0.0),VectorValue(0.0,0.0)]) +U = TrialFESpace(V,[uin,VectorValue(0.0,0.0),VectorValue(0.0,0.0),VectorValue(0.0,0.0)]) P = TrialFESpace(Q) R = TrialFESpace(T) @@ -90,20 +118,20 @@ Y = MultiFieldFESpace([V,Q,T]) # Properties Re = 60 # Reynolds number ρ = 1.0 # Density -L = 1.0 # Characteristic length +cl = H # Characteristic length u0_max = maximum(abs,get_dirichlet_dof_values(U)) -μ = ρ*L*u0_max/Re # Viscosity +μ = ρ*cl*u0_max/Re # Viscosity # Stabilization parameters γ = 1000.0 # Terms -σf_n(u,p) = μ*∇(u)⋅n_Γ - p*n_Γ +σf_n(u,p,n) = μ*∇(u) ⋅ n - p*n a_Ω(u,v) = μ*(∇(u) ⊙ ∇(v)) b_Ω(v,p) = - (∇⋅v)*p a_fluid((u,p),(v,q)) = - ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p)) * dΩ + - ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p) + (γ/h)*u⋅v ) * dΩout + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p)) * Ω.dΩf + + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p) + (γ/h)*u⋅v ) * Ω.dΩs ## Structure # Stabilization and material parameters @@ -112,30 +140,61 @@ function lame_parameters(E,ν) μ = E/(2*(1+ν)) (λ, μ) end -λs, μs = lame_parameters(1.0,0.3) +λs, μs = lame_parameters(0.1,0.05) #1.0,0.3) ϵ = (λs + 2μs)*1e-3 # Terms σ(ε) = λs*tr(ε)*one(ε) + 2*μs*ε -a_solid(d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))dΩout + - ∫(ϵ*(ε(s) ⊙ (σ ∘ ε(d))))dΩ # Ersatz +a_solid(d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))Ω.dΩs + + ∫(ϵ*(ε(s) ⊙ (σ ∘ ε(d))))Ω.dΩf # Ersatz ## Full problem -a((u,p,d),(v,q,s)) = a_fluid((u,p),(v,q)) + a_solid(d,s) + - ∫(σf_n(u,p) ⋅ s)dΓ # plus sign because of the normal direction -l((v,q,s)) = 0.0 - -op = AffineFEOperator(a,l,X,Y) - -uh, ph, dh = solve(op) - -# Mass flow rate through surface (this should be close to zero) -@show m = sum(∫(ρ*uh⋅n_Γ)dΓ) - -writevtk(Ω_act,path*"fsi-stokes-brinkmann-P2P1_elast-ersatz_full", - cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) -writevtk(Ω,path*"fsi-stokes-brinkmann-P2P1_elast-ersatz_fluid", +vec0 = VectorValue(0.0,0.0) +# minus sign because of the normal direction +function a_coupled((u,p,d),(v,q,s),φ) + n_AD = get_normal_vector(Ω.Γ) + return a_fluid((u,p),(v,q)) + a_solid(d,s) + + ∫(-σf_n(u,p,n_AD) ⋅ s)Ω.dΓ + ∫((-σf_n(u,p,vec0) ⋅ s))dΩ_act +end +l_coupled((v,q,s),φ) = ∫(0.0q)dΩ_act + +## Optimisation functionals +J_pres((u,p,d),φ) = ∫(p)dΓf_D - ∫(p)dΓf_N +J_comp((u,p,d),φ) = ∫(ε(d) ⊙ (σ ∘ ε(d)))Ω.dΩs +Vol((u,p,d),φ) = ∫(100*1/vol_D)Ω.dΩs - ∫(100*vf/vol_D)dΩ_act + +## Setup solver and FE operators +state_map = AffineFEStateMap(a_coupled,l_coupled,X,Y,V_φ,U_reg,φh) +pcfs = PDEConstrainedFunctionals(J_comp,[Vol],state_map) + +## Evolution Method +evo = CutFEMEvolve(V_φ,Ω,dΩ_act,h;max_steps) +reinit = StabilisedReinit(V_φ,Ω,dΩ_act,h;stabilisation_method=ArtificialViscosity(3.0)) +ls_evo = UnfittedFEEvolution(evo,reinit) + +## Hilbertian extension-regularisation problems +α = α_coeff*h_refine +a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ_act; +vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) + +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; + γ=γ_evo,verbose=true,constraint_names=[:Vol]) +for (it,(uh,ph,dh),φh) in optimiser + if iszero(it % iter_mod) + writevtk(Ω_act,path*"Omega_act_$it", + cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"ph"=>ph,"dh"=>dh]) + writevtk(Ω.Ωf,path*"Omega_f_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) + writevtk(Ω.Ωs,path*"Omega_s_$it", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) + writevtk(Ω.Γ,path*"Gamma_$it",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅Ω.n_Γ,"σf_n"=>σf_n(uh,ph,φh)]) + end + write_history(path*"/history.txt",optimiser.history) +end +it = get_history(optimiser).niter; uh,ph,dh = get_state(pcfs) +writevtk(Ω_act,path*"Omega_act_$it", + cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ω.Ωf,path*"Omega_f_$it", cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) -writevtk(Ωout,path*"fsi-stokes-brinkmann-P2P1_elast-ersatz_solid", +writevtk(Ω.Ωs,path*"Omega_s_$it", cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) - -writevtk(Γ,path*"fsi-stokes-brinkmann-P2P1_elast-ersatz_interface",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅n_Γ,"σf_n"=>σf_n(uh,ph)]) \ No newline at end of file +writevtk(Ω.Γ,path*"Gamma_$it",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅Ω.n_Γ,"σf_n"=>σf_n(uh,ph,φh)]) \ No newline at end of file diff --git a/src/Embedded/DifferentiableTriangulations.jl b/src/Embedded/DifferentiableTriangulations.jl index 3d40da7..9221a59 100644 --- a/src/Embedded/DifferentiableTriangulations.jl +++ b/src/Embedded/DifferentiableTriangulations.jl @@ -78,10 +78,10 @@ function Geometry.get_cell_reffe(t::DifferentiableTriangulation) get_cell_reffe(t.trian) end -# TODO: Do we ever need to dualize the cell points? +# TODO: Do we ever need to dualize the cell points? # I think its not necessary, since all the dual numbers are propagated through the cellmaps... -# Also: The current version dualizes only the phys points... -# If we want to indeed dualize this, we should probably also dualize the ref points +# Also: The current version dualizes only the phys points... +# If we want to indeed dualize this, we should probably also dualize the ref points # in the case where ttrian.trian is a SubCellTriangulation (but not in the case of a SubFacetTriangulation) # Anyway, I don't think this matters for now... function CellData.get_cell_points(ttrian::DifferentiableTriangulation) diff --git a/src/Embedded/IsolatedVolumes.jl b/src/Embedded/IsolatedVolumes.jl index 03e24ab..85a3c73 100644 --- a/src/Embedded/IsolatedVolumes.jl +++ b/src/Embedded/IsolatedVolumes.jl @@ -21,7 +21,7 @@ end ) Starting from a cell `cell`, crawls the cell graph provided by `cell_to_nbors` and colors all cells -connected to `cell` that +connected to `cell` that - belong to the group `group` (i.e., `cell_to_state[cell] ∈ group`), and - have not been seen yet (i.e., `!touched[cell]`). @@ -35,15 +35,15 @@ function tag_volume!( touched::BitVector ) where N @assert cell_to_state[cell] ∈ group - + q = Queue{Int}() enqueue!(q,cell) touched[cell] = true - + while !isempty(q) cell = dequeue!(q) cell_to_color[cell] = color - + nbors = cell_to_nbors[cell] for nbor in nbors state = cell_to_state[nbor] @@ -62,7 +62,7 @@ end groups = Tuple(unique(cell_to_state)) ) -Given a DiscreteModel `model` and an initial coloring `cell_to_state`, +Given a DiscreteModel `model` and an initial coloring `cell_to_state`, returns another coloring such that each color corresponds to a connected component of the graph of cells that are connected by a face and have their state in the same group. """ @@ -90,7 +90,7 @@ function tag_isolated_volumes( end cell += 1 end - + return cell_to_color, color_to_group end @@ -118,24 +118,31 @@ end """ get_isolated_volumes_mask(cutgeo::EmbeddedDiscretization,dirichlet_tags) -Given an EmbeddedDiscretization `cutgeo` and a list of tags `dirichlet_tags`, +Given an EmbeddedDiscretization `cutgeo` and a list of tags `dirichlet_tags`, this function returns a CellField which is `1` on isolated volumes and `0` otherwise. -We define an isolated volume as a volume that is IN but is not constrained by any +We define an isolated volume as a volume that is IN but is not constrained by any of the tags in `dirichlet_tags`. + +Specify the in domain using `IN_is`. Taking `IN_is = OUT` will find isolated +volumes for the `OUT` domain. """ function get_isolated_volumes_mask( - cutgeo::EmbeddedDiscretization,dirichlet_tags + cutgeo::EmbeddedDiscretization,dirichlet_tags; + IN_is = IN ) model = get_background_model(cutgeo) geo = get_geometry(cutgeo) bgcell_to_inoutcut = compute_bgcell_to_inoutcut(model,geo) - cell_to_color, color_to_group = tag_isolated_volumes(model,bgcell_to_inoutcut;groups=((CUT,IN),OUT)) + # Note: We switch IN and OUT based on the IN_is parameter. Namely, + # IN_is = IN -> groups = ((CUT,IN),IN*IN) = ((CUT,IN),OUT) + # IN_is = OUT -> groups = ((CUT,OUT),IN*OUT) = ((CUT,IN),IN) + cell_to_color, color_to_group = tag_isolated_volumes(model,bgcell_to_inoutcut;groups=((CUT,IN_is),IN*IN_is)) color_to_tagged = find_tagged_volumes(model,dirichlet_tags,cell_to_color,color_to_group) - + data = map(c -> !color_to_tagged[c] && isone(color_to_group[c]), cell_to_color) - return CellField(collect(Float64,data),Triangulation(model)) + return CellField(collect(Float64,data),Triangulation(model)) end # Distributed @@ -158,13 +165,14 @@ function tag_isolated_volumes( end function get_isolated_volumes_mask( - cutgeo::GridapEmbedded.Distributed.DistributedEmbeddedDiscretization,dirichlet_tags + cutgeo::GridapEmbedded.Distributed.DistributedEmbeddedDiscretization,dirichlet_tags; + IN_is = IN ) model = get_background_model(cutgeo) geo = get_geometry(cutgeo) bgcell_to_inoutcut = map(compute_bgcell_to_inoutcut,local_views(model),local_views(geo)) - cell_to_lcolor, lcolor_to_group, color_gids = tag_isolated_volumes(model,bgcell_to_inoutcut,((CUT,IN),OUT)) + cell_to_lcolor, lcolor_to_group, color_gids = tag_isolated_volumes(model,bgcell_to_inoutcut,((CUT,IN_is),IN*IN_is)) lcolor_to_tagged = map(local_views(model),cell_to_lcolor,lcolor_to_group) do model, cell_to_lcolor, lcolor_to_group find_tagged_volumes(model,dirichlet_tags,cell_to_lcolor,lcolor_to_group) @@ -176,7 +184,7 @@ function get_isolated_volumes_mask( trian = Triangulation(GridapDistributed.WithGhost(),model) fields = map(local_views(trian),cell_to_lcolor,lcolor_to_group,lcolor_to_tagged) do trian, cell_to_lcolor, lcolor_to_group, lcolor_to_tagged data = map(c -> !lcolor_to_tagged[c] && isone(lcolor_to_group[c]), cell_to_lcolor) - CellField(collect(Float64,data),trian) + CellField(collect(Float64,data),trian) end return GridapDistributed.DistributedCellField(fields,trian) end @@ -192,7 +200,7 @@ function generate_volume_gids( cache.buffer_snd[k] = cell_to_lcolor[lid] end cache.neighbors_snd,cache.neighbors_rcv,cache.buffer_snd,cache.buffer_rcv - end |> tuple_of_arrays; + end |> tuple_of_arrays; graph = ExchangeGraph(neighbors_snd,neighbors_rcv) t = exchange!(buffer_rcv,buffer_snd,graph) diff --git a/src/GridapExtensions.jl b/src/GridapExtensions.jl index 45cdea0..351199c 100644 --- a/src/GridapExtensions.jl +++ b/src/GridapExtensions.jl @@ -32,7 +32,7 @@ end # Transpose contributions before assembly -transpose_contributions(b::DistributedDomainContribution) = +transpose_contributions(b::DistributedDomainContribution) = DistributedDomainContribution(map(transpose_contributions,local_views(b))) function transpose_contributions(b::DomainContribution) diff --git a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl index b72b20d..8cda263 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl +++ b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl @@ -53,20 +53,22 @@ function get_transient_operator(φh,velh,s::CutFEMEvolve) γg, h, dΓg, n_Γg = params.γg, params.h, params.dΓg, params.n_Γg ϵ = 1e-20 - β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh)) + v_norm = maximum(abs,get_free_dof_values(velh)) + β = velh/(ϵ + v_norm) * ∇(φh)/(ϵ + norm ∘ ∇(φh)) + stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ_bg + ∫(γg*h^2*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ_bg forcing(t,v) = ∫(0v)dΩ_bg + ∫(0*jump(∇(v) ⋅ n_Γg))dΓg - # Second term is added to address the following issue: + # Second term is added to address the following issue: # - ODEs is allocating separately the residual and jacobian - # - This is fine in serial, but in parallel there are some instances where the the following happens: + # - This is fine in serial, but in parallel there are some instances where the the following happens: # - The residual is touched by less ghost entries than the columns of the matrix - # - If we assemble both jac and res together, we communicate the extra ghost ids to - # the residual, so everything is consistent. - # - However, if we assemble the residual and jacobian separately, + # - If we assemble both jac and res together, we communicate the extra ghost ids to + # the residual, so everything is consistent. + # - However, if we assemble the residual and jacobian separately, # the residual is not aware of the extra ghost ids - # This happens when there are touched ghost entries that do not belong to the local domain. - # In particular, this happens when we have jumps, where some contributions come from two + # This happens when there are touched ghost entries that do not belong to the local domain. + # In particular, this happens when we have jumps, where some contributions come from two # cells away. Boundary cells then get contributions from cells which are not in the local domain. Ut_φ = TransientTrialFESpace(V_φ) ode_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ; @@ -91,10 +93,9 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing) h, max_steps = params.h, params.max_steps # Setup FE operator and solver + ode_solver.dt = γ*h ode_op = get_transient_operator(φh,velh,s) - 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) + ode_sol = solve(ode_solver,ode_op,0.0,ode_solver.dt*max_steps,φh) # March march = Base.iterate(ode_sol) @@ -124,6 +125,7 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache) # We do this so that the first iterate of ODESolution always recomputes the # stiffness matrix and associated the Jacboian, numerical setups, etc via # `constant_forms = (false,true)`. + ode_solver.dt = γ*h ode_op = get_transient_operator(φh,velh,s) # Between the first iterate and subsequent iterates we use the function # `update_reuse!` to update the iterator state so that we re-use @@ -131,11 +133,9 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache) # whether we are solving a new ODE with the same functional form but # updated coefficients in the weak form. If so, we want to re-use the cache. state_inter = update_reuse!(cache,false;zero_tF=true) - 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) + ode_sol = solve(ode_solver,ode_op,0.0,ode_solver.dt*max_steps,φh) march = Base.iterate(ode_sol,state_inter) # First step includes stiffness matrix update data, state = march state_updated = update_reuse!(state,true) # Fix the stiffness matrix for remaining march @@ -151,10 +151,4 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache) s.cache = state_updated update_collection!(s.Ωs,φh) # TODO: remove? 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/CutFEMEvolve_new.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve_new.jl index d1cb215..649bd13 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve_new.jl +++ b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve_new.jl @@ -18,46 +18,42 @@ Burman et al. (2017). DOI: `10.1016/j.cma.2017.09.005`. - The stepsize `dt = 0.1` in `RungeKutta` is a place-holder and is updated using the `γ` passed to `solve!`. """ -mutable struct CutFEMEvolve{A,B} <: Evolver +mutable struct CutFEMEvolve{A,B,C} <: Evolver ode_solver::ODESolver Ωs::EmbeddedCollection - space::A + dΩ_bg::A + space::B assembler::Assembler - params::B + params::C cache - function CutFEMEvolve(V_φ::A,Ωs::EmbeddedCollection,h::Real; + function CutFEMEvolve(V_φ::B,Ωs::EmbeddedCollection,dΩ_bg::A,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_φ), - Γg_quad_degree = 2get_order(V_φ), - Ωφ_quad_degree = 2get_order(V_φ)) where A - + assembler=SparseMatrixAssembler(V_φ,V_φ)) where {A,B} Γg = SkeletonTriangulation(get_triangulation(V_φ)) - dΓg = Measure(Γg,Γg_quad_degree) + dΓg = Measure(Γg,2get_order(V_φ)) n_Γg = get_normal_vector(Γg) - Ωφ = get_triangulation(V_φ) - dΩφ = Measure(Ωφ,Ωφ_quad_degree) - params = (;γg,h,max_steps,dΓg,n_Γg,dΩφ) - new{A,typeof(params)}(ode_solver,Ωs,V_φ,assembler,params,nothing) + params = (;γg,h,max_steps,dΓg,n_Γg) + new{A,B,typeof(params)}(ode_solver,Ωs,dΩ_bg,V_φ,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_measure(s::CutFEMEvolve) = get_params(s).dΩφ +get_measure(s::CutFEMEvolve) = s.dΩ_bg get_params(s::CutFEMEvolve) = s.params get_cache(s::CutFEMEvolve) = s.cache -function get_transient_operator(φh::CellField,velh::CellField,s::CutFEMEvolve) - V_φ, assembler, params = s.space, s.assembler, s.params - γg, h, dΓg, n_Γg, dΩ_bg = params.γg, params.h, params.dΓg, params.n_Γg, params.dΩφ +function get_transient_operator(φh,velh,s::CutFEMEvolve,dt) + V_φ, dΩ_bg, assembler, params = s.space, s.dΩ_bg, s.assembler, s.params + γg, h, dΓg, n_Γg = params.γg, params.h, params.dΓg, params.n_Γg ϵ = 1e-20 - β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh)) + β = dt*velh*∇(φh)/(ϵ + norm ∘ ∇(φh)) stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ_bg + ∫(γg*h^2*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ_bg forcing(t,v) = ∫(0v)dΩ_bg + ∫(0*jump(∇(v) ⋅ n_Γg))dΓg @@ -89,16 +85,16 @@ function update_reuse!(state,reuse_new;zero_tF=false) return U, (_tF, stateF, state0, uF, odecache_new) end -function solve!(s::CutFEMEvolve,φh::CellField,velh::CellField,γ,cache::Nothing) +function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing) ode_solver = get_ode_solver(s) params = get_params(s) h, max_steps = params.h, params.max_steps # Setup FE operator and solver - ode_op = get_transient_operator(φh,velh,s) 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,interpolate(φh,s.space)) + ode_solver.dt = γ*h + ode_op = get_transient_operator(φh,velh,s,dt/(γ*h)) + ode_sol = solve(ode_solver,ode_op,0.0,ode_solver.dt*max_steps,φh) # March march = Base.iterate(ode_sol) @@ -112,19 +108,13 @@ function solve!(s::CutFEMEvolve,φh::CellField,velh::CellField,γ,cache::Nothing # Update φh and cache _, φhF = data - if get_triangulation(φh) === get_triangulation(φhF) - copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) - else - interpolate!(φhF,get_free_dof_values(φh),get_fe_space(φh)) - end - writevtk(get_triangulation(φhF),"./results/FCM_thermal_compliance_ALM/test",cellfields=["φ"=>φhF]) - writevtk(get_triangulation(φh),"./results/FCM_thermal_compliance_ALM/test2",cellfields=["φ"=>φh]) + copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) s.cache = state_new update_collection!(s.Ωs,φh) # TODO: remove? return φh end -function solve!(s::CutFEMEvolve,φh::CellField,velh::CellField,γ,cache) +function solve!(s::CutFEMEvolve,φh,velh,γ,cache) ode_solver = get_ode_solver(s) params = get_params(s) h, max_steps = params.h, params.max_steps @@ -134,18 +124,18 @@ function solve!(s::CutFEMEvolve,φh::CellField,velh::CellField,γ,cache) # We do this so that the first iterate of ODESolution always recomputes the # stiffness matrix and associated the Jacboian, numerical setups, etc via # `constant_forms = (false,true)`. - ode_op = get_transient_operator(φh,velh,s) + dt = _compute_Δt(h,γ,get_free_dof_values(velh)) + ode_solver.dt = γ*h + ode_op = get_transient_operator(φh,velh,s,dt/(γ*h)) # Between the first iterate and subsequent iterates we use the function # `update_reuse!` to update the iterator state so that we re-use # the stiffness matrix, etc. The Optional argument `zero_tF` indicates # whether we are solving a new ODE with the same functional form but # updated coefficients in the weak form. If so, we want to re-use the cache. state_inter = update_reuse!(cache,false;zero_tF=true) - 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,interpolate(φh,s.space)) + ode_sol = solve(ode_solver,ode_op,0.0,ode_solver.dt*max_steps,φh) march = Base.iterate(ode_sol,state_inter) # First step includes stiffness matrix update data, state = march state_updated = update_reuse!(state,true) # Fix the stiffness matrix for remaining march @@ -157,11 +147,7 @@ function solve!(s::CutFEMEvolve,φh::CellField,velh::CellField,γ,cache) ## Update φh and cache _, φhF = data - if get_triangulation(φh) === get_triangulation(φhF) - copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) - else - interpolate!(φhF,get_free_dof_values(φh),get_fe_space(φh)) - end + copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) s.cache = state_updated update_collection!(s.Ωs,φh) # TODO: remove? return φh diff --git a/test/seq/EmbeddedDifferentiationTests.jl b/test/seq/EmbeddedDifferentiationTests.jl index cd12b39..f84d95d 100644 --- a/test/seq/EmbeddedDifferentiationTests.jl +++ b/test/seq/EmbeddedDifferentiationTests.jl @@ -26,6 +26,8 @@ function level_set(shape::Symbol;N=4) x -> abs(x[1]-0.5)+abs(x[2]-0.5)-0.25-0/n/10 # Diamond elseif shape == :circle x -> sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.5223 # Circle + elseif shape == :circle_2 + x -> sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.23 # Circle elseif shape == :square_prism x -> max(abs(x[1]-0.5),abs(x[2]-0.5),abs(x[3]-0.5))-0.25 # Square prism elseif shape == :corner_3d @@ -198,20 +200,97 @@ function main( end end +## Concering integrals of the form `φ->∫(f ⋅ n(φ))dΓ(φ)` +function main_normal( + model,φ::Function,f::Function; + vtk=false, + name="embedded", + run_test=true +) + order = 1 + reffe = ReferenceFE(lagrangian,Float64,order) + V_φ = TestFESpace(model,reffe) + + φh = interpolate(φ,V_φ) + + geo = DiscreteGeometry(φh,model) + cutgeo = cut(model,geo) + + Γ = EmbeddedBoundary(cutgeo) + Γ_AD = DifferentiableTriangulation(Γ,V_φ) + dΓ_AD = Measure(Γ_AD,2*order) + dΓ = Measure(Γ,2*order) + + fh_Γ = CellField(f,Γ) + fh_Γ_AD = CellField(f,Γ_AD) + + function J_int(φ) + n = get_normal_vector(Γ_AD) + ∫(fh_Γ_AD⋅n)dΓ_AD + end + dJ_int_AD = gradient(J_int,φh) + dJ_int_AD_vec = assemble_vector(dJ_int_AD,V_φ) + + _n(∇φ) = ∇φ/(10^-20+norm(∇φ)) + dJ_int_phi = ∇(φ->∫(fh_Γ_AD ⋅ (_n ∘ (∇(φ))))dΓ_AD,φh) + dJh_int_phi = assemble_vector(dJ_int_phi,V_φ) + + run_test && @test norm(dJ_int_AD_vec - dJh_int_phi) < 1e-10 + + # Analytic + # Note: currently, the analytic result is only valid on closed domains thanks + # to the divergence theorem. I think it would take significant work to compute + # the analytic derivative generally as we can't rely on divergence theorem to + # rewrite it in a convenient way. As a result, we don't have an analytic result + # for general cases such as ∫( f(n(φ)) )dΓ(φ), nor the case when Γ intersects + # ∂D. Thankfully, we have AD instead ;) + # Note 2: For the case that Γ does intersect the surface, the result is correct + # everywhere except on the intersection. + + fh2(x) = VectorValue((1-x[1])^2,(1-x[2])^2) + fh_Γ = CellField(fh2,Γ) + fh_Γ_AD = CellField(fh2,Γ_AD) + + # Note: this comes from rewriting via the divergence theorem: + # ∫(f ⋅ n(φ))dΓ(φ) = ∫(∇⋅f)dΩ(φ) + dJ_int_exact3(w) = ∫(-(∇⋅(fh_Γ))*w/(norm ∘ (∇(φh))))dΓ + dJh_int_exact3 = assemble_vector(dJ_int_exact3,V_φ) + + run_test && @test norm(dJh_int_exact3 - dJ_int_AD_vec) < 1e-10 + + if vtk + path = "results/$(name)" + Ω_bg = Triangulation(model) + writevtk(Ω_bg,path,cellfields=[ + "dJ_AD"=>FEFunction(V_φ,dJ_int_AD_vec), + "dJ_AD_with_phi"=>FEFunction(V_φ,dJh_int_phi), + "dJ_exact"=>FEFunction(V_φ,dJh_int_exact3) + ]) + end +end + ####################### D = 2 n = 10 model = generate_model(D,n) φ = level_set(:circle) -f = x -> 1.0 +f(x) = 1.0 main(model,φ,f;vtk=true) D = 2 n = 10 model = generate_model(D,n) φ = level_set(:circle) -f = x -> x[1]+x[2] +f(x) = x[1]+x[2] main(model,φ,f;vtk=true) +D = 2 +n = 10 +model = generate_model(D,n) +φ = level_set(:circle_2) +fvec((x,y)) = VectorValue((1-x)^2,(1-y)^2) +main_normal(model,φ,fvec;vtk=true) +# main_normal(model,level_set(:circle),fvec;vtk=true) # This will fail as expected + end \ No newline at end of file