Add gmsh test (WIP)

Project.toml

@@ -15,6 +15,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapDistributed = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
 GridapEmbedded = "8838a6a3-0006-4405-b874-385995508d5d"
+GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
 GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
@@ -33,6 +34,7 @@ DataStructures = "0.18.20"
 FillArrays = "1"
 Gridap = "0.18"
 GridapDistributed = "0.4"
+GridapGmsh = "0.7.2"
 GridapPETSc = "0.5"
 GridapSolvers = "0.4"
 JLD2 = "0.4,0.5"

scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl

@@ -0,0 +1,189 @@
+using Gridap, Gridap.Geometry, Gridap.Adaptivity
+using GridapEmbedded, GridapEmbedded.LevelSetCutters
+using GridapGmsh
+using GridapTopOpt
+
+path = "./results/GMSH-TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi/"
+mkpath(path)
+
+γ_evo = 0.1
+max_steps = 20
+vf = 0.03
+α_coeff = 2
+iter_mod = 1
+D = 2
+
+# Load gmsh mesh (Currently need to update mesh.geo and these values concurrently)
+L = 2.0;
+H = 0.5;
+x0 = 0.5;
+l = 0.4;
+w = 0.05;
+a = 0.3;
+b = 0.03;
+
+model = GmshDiscreteModel((@__DIR__)*"/mesh.msh")
+# writevtk(model,path*"model")
+
+# el_Δ = get_el_Δ(_model)
+# h = maximum(el_Δ)
+# h_refine = maximum(el_Δ)/2
+
+# Cut the background model
+reffe_scalar = ReferenceFE(lagrangian,Float64,1)
+V_φ = TestFESpace(model,reffe_scalar)
+V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Omega_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
+
+Ω_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
+
+writevtk(Ω_act,path*"Omega_act",
+  cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh))])
+writevtk(Ω.Ωf,path*"Omega_f")
+writevtk(Ω.Ωs,path*"Omega_s")
+writevtk(Ω.Γ,path*"Gamma")
+
+# 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-1)
+
+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_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))
+μ = ρ*cl*u0_max/Re # Viscosity
+# Stabilization parameters
+γ = 1000.0
+
+# Terms
+σ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Ωf +
+  ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p) + (γ/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(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Ω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_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(Ω.Ω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)])

scripts/Embedded/Examples/fsi/gmsh/mesh.geo

@@ -0,0 +1,111 @@
+// Gmsh project created on Fri Dec 06 15:33:11 2024
+//+
+SetFactory("OpenCASCADE");
+Mesh.SaveAll=1;
+
+//+ 8-----7---------6-----------5
+//+ |     |         |           |
+//+ |     9--13-14--10		      |
+//+ |     |  |   |  |           |
+//+ |     11-12 15-16           |
+//+ 1-----2---------3-----------4
+//+ Area enclosed by 2,3,16,15,14,13,12,11,2 is nondesign domain
+
+//+ Dims
+L = 2.0;
+H = 0.5;
+x0 = 0.5;
+l = 0.4;
+w = 0.05;
+a = 0.3;
+b = 0.03;
+
+//+ Mesh sizes
+size_f = 0.05;
+size_s = 0.005;
+
+//+ Main area
+Point(1) = {0,0,0,size_f};
+Point(2) = {x0-l/2,0,0,size_s};
+Point(3) = {x0+l/2,0,0,size_s};
+Point(4) = {L,0,0,size_f};
+Point(5) = {L,H,0,size_f};
+Point(6) = {x0+l/2,H,0,size_f};
+Point(7) = {x0-l/2,H,0,size_f};
+Point(8) = {0,H,0,size_f};
+Point(9) = {x0-l/2,a,0,size_s};
+Point(10) = {x0+l/2,a,0,size_s};
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,5};
+Line(5) = {5,6};
+Line(6) = {6,7};
+Line(7) = {7,8};
+Line(8) = {8,1};
+
+//+ Non-designable region
+Point(11) = {x0-l/2,b,0,size_s};
+Point(12) = {x0-w/2,b,0,size_s};
+Point(13) = {x0-w/2,a,0,size_s};
+Point(14) = {x0+w/2,a,0,size_s};
+Point(15) = {x0+w/2,b,0,size_s};
+Point(16) = {x0+l/2,b,0,size_s};
+
+Line(9) = {7,9};
+Line(10) = {6,10};
+Line(11) = {9,11};
+Line(12) = {10,16};
+Line(13) = {9,13};
+Line(14) = {14,10};
+
+Line(15) = {2,11};
+Line(16) = {11,12};
+Line(17) = {12,13};
+Line(18) = {13,14};
+Line(19) = {14,15};
+Line(20) = {15,16};
+Line(21) = {16,3};
+
+//+
+Curve Loop(1) = {7, 8, 1, 15, -11, -9};
+//+
+Plane Surface(1) = {1};
+//+
+Curve Loop(2) = {6, 9, 13, 18, 14, -10};
+//+
+Plane Surface(2) = {2};
+//+
+Curve Loop(3) = {5, 10, 12, 21, 3, 4};
+//+
+Plane Surface(3) = {3};
+//+
+Curve Loop(4) = {13, -17, -16, -11};
+//+
+Plane Surface(4) = {4};
+//+
+Curve Loop(5) = {19, 20, -12, -14};
+//+
+Plane Surface(5) = {5};
+//+
+Curve Loop(6) = {18, 19, 20, 21, -2, 15, 16, 17};
+//+ Non-design domain
+Plane Surface(6) = {6};
+
+//+
+Physical Curve("Gamma_f_D", 22) = {8};
+Physical Point("Gamma_f_D", 23) = {8, 1};
+//+
+Physical Curve("Gamma_f_N", 24) = {4};
+Physical Point("Gamma_f_N", 25) = {5, 4};
+//+
+Physical Curve("Gamma_NoSlipTop", 26) = {7, 6, 5};
+Physical Point("Gamma_NoSlipTop", 27) = {7, 6};
+Physical Curve("Gamma_NoSlipBottom", 28) = {1, 3};
+//+
+Physical Curve("Gamma_s_D", 29) = {2};
+Physical Point("Gamma_s_D", 30) = {2, 3};
+//+
+Physical Surface("Omega_NonDesign", 31) = {6};
+Physical Curve("Omega_NonDesign", 32) = {17, 18, 19, 20, 16, 15, 21};