From e5d0fce95af0c5b376d86830a51aee6b150c0620 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Sun, 1 Dec 2024 16:23:12 +1000 Subject: [PATCH] Updated MWEs --- .../CutFEM_2d_thermal_(island_checking).jl | 120 ++++++++++++++++ .../fsi/2-cutfem_navier-stokes_fsi.jl | 2 +- .../Examples/fsi/2-cutfem_stokes_fsi.jl | 2 +- ...Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl | 132 ++++++++++++++++++ scripts/Embedded/MWEs/isolated_volumes.jl | 22 +-- test/seq/IsolatedVolumeTests.jl | 1 + 6 files changed, 267 insertions(+), 12 deletions(-) create mode 100644 scripts/Embedded/Examples/CutFEM_2d_thermal_(island_checking).jl create mode 100644 scripts/Embedded/Examples/fsi/5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl create mode 100644 test/seq/IsolatedVolumeTests.jl diff --git a/scripts/Embedded/Examples/CutFEM_2d_thermal_(island_checking).jl b/scripts/Embedded/Examples/CutFEM_2d_thermal_(island_checking).jl new file mode 100644 index 0000000..cb18fa6 --- /dev/null +++ b/scripts/Embedded/Examples/CutFEM_2d_thermal_(island_checking).jl @@ -0,0 +1,120 @@ +using Gridap,GridapTopOpt, GridapSolvers +using Gridap.Adaptivity, Gridap.Geometry +using GridapEmbedded, GridapEmbedded.LevelSetCutters + +using GridapTopOpt: StateParamIntegrandWithMeasure + +path="./results/CutFEM_thermal_compliance_ALM_island_detect/" +rm(path,force=true,recursive=true) +mkpath(path) +n = 50 +order = 1 +γ = 0.1 +max_steps = floor(Int,order*n/5) +vf = 0.4 +α_coeff = 4max_steps*γ +iter_mod = 1 + +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +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.0 && (x[2] <= 0.2 + eps() || x[2] >= 0.8 - eps())) +f_Γ_N(x) = (x[1] ≈ 1 && 0.4 - eps() <= x[2] <= 0.6 + eps()) +update_labels!(1,model,f_Γ_D,"Gamma_D") +update_labels!(2,model,f_Γ_N,"Gamma_N") + +## Triangulations and measures +Ω = Triangulation(model) +Γ_N = BoundaryTriangulation(model,tags="Gamma_N") +dΩ = Measure(Ω,2*order) +dΓ_N = Measure(Γ_N,2*order) +vol_D = sum(∫(1)dΩ) + +## Levet-set function space and derivative regularisation space +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"]) +U_reg = TrialFESpace(V_reg,0) +V_φ = TestFESpace(model,reffe_scalar) +V_χ = TestFESpace(model,ReferenceFE(lagrangian,Float64,0)) + +## Levet-set function +φh = interpolate(x->-cos(4π*x[1])*cos(4π*x[2])-0.4,V_φ) +Ωs = EmbeddedCollection(model,φh) do cutgeo,_ + Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL),V_φ) + Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ) + Γg = GhostSkeleton(cutgeo) + Ωact = Triangulation(cutgeo,ACTIVE) + (; + :Ωin => Ωin, + :dΩin => Measure(Ωin,2*order), + :Γg => Γg, + :dΓg => Measure(Γg,2*order), + :n_Γg => get_normal_vector(Γg), + :Γ => Γ, + :dΓ => Measure(Γ,2*order), + :Ωact => Ωact, + :χ => GridapTopOpt.get_isolated_volumes_mask(cutgeo,["Gamma_D"]) + ) +end + +## Weak form +const γg = 0.1 +a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + + ∫((γg*h)*jump(Ωs.n_Γg⋅∇(v))*jump(Ωs.n_Γg⋅∇(u)))Ωs.dΓg + + ∫(Ωs.χ*v*u)Ωs.dΩin +l(v,φ) = ∫(v)dΓ_N + +## Optimisation functionals +J(u,φ) = a(u,u,φ) +Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ +dVol(q,u,φ) = ∫(-1/vol_D*q/(1e-20 + norm ∘ (∇(φ))))Ωs.dΓ + +## Setup solver and FE operators +state_collection = EmbeddedCollection(model,φh) do _,_ + V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,0.0) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) + (; + :state_map => state_map, + :J => StateParamIntegrandWithMeasure(J,state_map), + :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,]) + ) +end +pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,)) + +## Evolution Method +evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps) +reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0))#InteriorPenalty(V_φ)) +ls_evo = UnfittedFEEvolution(evo,reinit) +reinit!(ls_evo,φh) + +## Hilbertian extension-regularisation problems +α = α_coeff*(h_refine/order)^2 +a_hilb(p,q) =∫(α*∇(p)⋅∇(q) + p*q)dΩ; +vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) + +## Optimiser +converged(m) = GridapTopOpt.default_al_converged( + m; + L_tol = 0.01*h_refine, + C_tol = 0.01 +) +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;debug=true, + γ,verbose=true,constraint_names=[:Vol],converged) +for (it,uh,φh,state) in optimiser + x_φ = get_free_dof_values(φh) + idx = findall(isapprox(0.0;atol=10^-10),x_φ) + !isempty(idx) && @warn "Boundary intersects nodes!" + if iszero(it % iter_mod) + writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"velh"=>FEFunction(V_φ,state.vel),"χ"=>Ωs.χ]) + writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh]) + end + write_history(path*"/history.txt",optimiser.history) +end +it = get_history(optimiser).niter; uh = get_state(pcfs) +writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"χ"=>Ωs.χ]) +writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh]) \ No newline at end of file diff --git a/scripts/Embedded/Examples/fsi/2-cutfem_navier-stokes_fsi.jl b/scripts/Embedded/Examples/fsi/2-cutfem_navier-stokes_fsi.jl index f3439de..8b93c0b 100644 --- a/scripts/Embedded/Examples/fsi/2-cutfem_navier-stokes_fsi.jl +++ b/scripts/Embedded/Examples/fsi/2-cutfem_navier-stokes_fsi.jl @@ -60,7 +60,7 @@ dΓi = Measure(Γi,degree) # Setup FESpace -uin(x) = VectorValue(x[2],0.0) +uin(x) = VectorValue(x[2]*(1-x[2]),0.0) reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P) reffe_p = ReferenceFE(lagrangian,Float64,order,space=:P) diff --git a/scripts/Embedded/Examples/fsi/2-cutfem_stokes_fsi.jl b/scripts/Embedded/Examples/fsi/2-cutfem_stokes_fsi.jl index 8e7e6b7..914ae17 100644 --- a/scripts/Embedded/Examples/fsi/2-cutfem_stokes_fsi.jl +++ b/scripts/Embedded/Examples/fsi/2-cutfem_stokes_fsi.jl @@ -58,7 +58,7 @@ dΓi = Measure(Γi,degree) # Setup FESpace -uin(x) = VectorValue(x[2],0.0) +uin(x) = VectorValue(x[2]*(1-x[2]),0.0) reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P) reffe_p = ReferenceFE(lagrangian,Float64,order,space=:P) diff --git a/scripts/Embedded/Examples/fsi/5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl b/scripts/Embedded/Examples/fsi/5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl new file mode 100644 index 0000000..4b1cb1c --- /dev/null +++ b/scripts/Embedded/Examples/fsi/5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl @@ -0,0 +1,132 @@ +using Gridap, Gridap.Geometry, Gridap.Adaptivity +using GridapEmbedded, GridapEmbedded.LevelSetCutters +using GridapTopOpt + +path = "./results/fsi testing/" +mkpath(path) + +# Cut the background model +n = 100 +partition = (n,n) +D = length(partition) +_model = CartesianDiscreteModel((0,1,0,1),partition) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model + +el_Δ = get_el_Δ(_model) +h = maximum(el_Δ) +f_Γ_D(x) = x[1] ≈ 0 +f_Γ_NoSlipTop(x) = x[2] ≈ 1 +f_Γ_NoSlipBottom(x) = x[2] ≈ 0 +update_labels!(1,model,f_Γ_D,"Gamma_D") +update_labels!(2,model,f_Γ_NoSlipTop,"Gamma_NoSlipTop") +update_labels!(3,model,f_Γ_NoSlipBottom,"Gamma_NoSlipBottom") + +# 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 +order = 1 +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 + +uin(x) = VectorValue(x[2]*(1-x[2]),0.0) # Change this!! + +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_D","Gamma_NoSlipTop","Gamma_NoSlipBottom"]) +Q = TestFESpace(Ω_act,reffe_p,conformity=:H1) +T = TestFESpace(Ω_act ,reffe_d,conformity=:H1,dirichlet_tags=["Gamma_NoSlipBottom"]) + +U = TrialFESpace(V,[uin,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 +L = 1.0 # Characteristic length +u0_max = maximum(abs,get_dirichlet_dof_values(U)) +μ = ρ*L*u0_max/Re # Viscosity +# Stabilization parameters +β1 = 0.2 +γ = 1000.0 + +# Terms +σf_n(u,p) = μ*∇(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Ω + + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p)-c_Ω(p,q) + (γ/h)*u⋅v ) * dΩout + +## 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) +ϵ = (λs + 2μs)*1e-3 +# Terms +σ(ε) = λs*tr(ε)*one(ε) + 2*μs*ε +a_solid(d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))dΩout + + ∫(ϵ*(ε(s) ⊙ (σ ∘ ε(d))))dΩ # 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_elast-ersatz_full", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ω,path*"fsi-stokes-brinkmann_elast-ersatz_fluid", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ωout,path*"fsi-stokes-brinkmann_elast-ersatz_solid", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) + +writevtk(Γ,path*"fsi-stokes-brinkmann_elast-ersatz_interface",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅n_Γ,"σf_n"=>σf_n(uh,ph)]) \ No newline at end of file diff --git a/scripts/Embedded/MWEs/isolated_volumes.jl b/scripts/Embedded/MWEs/isolated_volumes.jl index d875391..1efd8a0 100644 --- a/scripts/Embedded/MWEs/isolated_volumes.jl +++ b/scripts/Embedded/MWEs/isolated_volumes.jl @@ -6,20 +6,22 @@ using GridapEmbedded.LevelSetCutters using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Adaptivity, Gridap.Arrays order = 1 -n = 100 +n = 41 model = UnstructuredDiscreteModel(CartesianDiscreteModel((0,1,0,1),(n,n))) +update_labels!(1,model,x->x[1]≈0,"Gamma_D") Ω = Triangulation(model) reffe_scalar = ReferenceFE(lagrangian,Float64,order) V_φ = TestFESpace(model,reffe_scalar) -φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])-0.11,V_φ) - -# R = 0.195 -# R = 0.2 # This fails -# f(x0,r) = x -> sqrt((x[1]-x0[1])^2 + (x[2]-x0[2])^2) - r -# φh = interpolate(x->-f([0.5,0.5],R)(x),V_φ) -# φh = interpolate(x->min(f([0.25,0.5],R)(x),f([0.75,0.5],R)(x)),V_φ) +f(x,y0) = abs(x[2]-y0) - 0.05 +g(x,x0,y0,r) = sqrt((x[1]-x0)^2 + (x[2]-y0)^2) - r +f(x) = min(f(x,0.75),f(x,0.25), + g(x,0.15,0.5,0.1), + g(x,0.5,0.6,0.2), + g(x,0.85,0.5,0.1), + g(x,0.5,0.15,0.05)) +φh = interpolate(f,V_φ) geo = DiscreteGeometry(φh,model) cutgeo = cut(model,geo) @@ -27,10 +29,10 @@ cutgeo = cut(model,geo) bgcell_to_inoutcut = compute_bgcell_to_inoutcut(model,geo) cell_to_color, color_to_group = GridapTopOpt.tag_isolated_volumes(model,bgcell_to_inoutcut;groups=((GridapTopOpt.CUT,IN),OUT)) -color_to_tagged = GridapTopOpt.find_tagged_volumes(model,["tag_5","tag_7"],cell_to_color,color_to_group) +color_to_tagged = GridapTopOpt.find_tagged_volumes(model,["Gamma_D"],cell_to_color,color_to_group) cell_to_tagged = map(c -> color_to_tagged[c], cell_to_color) -μ = GridapTopOpt.get_isolated_volumes_mask(cutgeo,["tag_5","tag_7"]) +μ = GridapTopOpt.get_isolated_volumes_mask(cutgeo,["Gamma_D"]) writevtk( Ω,"results/background", diff --git a/test/seq/IsolatedVolumeTests.jl b/test/seq/IsolatedVolumeTests.jl new file mode 100644 index 0000000..f87f5c1 --- /dev/null +++ b/test/seq/IsolatedVolumeTests.jl @@ -0,0 +1 @@ +# TODO \ No newline at end of file