Merge branch 'feature-gridap-embedded-compat' of github.com:zjwegert/…

…GridapTopOpt.jl into feature-gridap-embedded-compat

scripts/Embedded/Examples/FCM_2d_thermal_with_island_checking.jl

@@ -8,7 +8,7 @@ using DataStructures
 
 const CUT = 0
 
-# TODO: Can be optimized CartesianModels
+# TODO: Can be optimized for CartesianModels
 function generate_neighbor_graph(model::DiscreteModel{Dc}) where Dc
   topo = get_grid_topology(model)
   cell_to_node = Geometry.get_faces(topo, Dc, 0)
@@ -19,22 +19,27 @@ function generate_neighbor_graph(model::DiscreteModel{Dc}) where Dc
   return cell_to_nbors
 end
 
-"""
-  Given an initial interface cell, enqueue all the CUT cells in the same interface
-  inside the queue `q_cut` and mark them as touched in the `touched` array.
-"""
-function enqueue_interface!(q_cut,cell_to_nbors,cell_to_inoutcut,touched,cell)
-  q = Queue{Int}(); enqueue!(q,cell)
-  enqueue!(q_cut,cell)
+function tag_volume!(
+  cell::Int,color::Int16,
+  cell_to_nbors::Vector{Vector{Int32}},
+  cell_to_inoutcut::Vector{Int8},
+  cell_to_color::Vector{Int16},
+  touched::BitVector
+)
+  q = Queue{Int}()
+  enqueue!(q,cell)
   touched[cell] = true
+  state = cell_to_inoutcut[cell]
+
   while !isempty(q)
     cell = dequeue!(q)
+    cell_to_color[cell] = color
+
     nbors = cell_to_nbors[cell]
     for nbor in nbors
-      if !touched[nbor] && (cell_to_inoutcut[nbor] == CUT)
-        touched[nbor] = true
-        enqueue!(q_cut,nbor)
+      if !touched[nbor] && (cell_to_inoutcut[nbor] == state)
         enqueue!(q,nbor)
+        touched[nbor] = true
       end
     end
   end
@@ -47,62 +52,42 @@ function tag_isolated_volumes(
   n_cells = num_cells(model)
   cell_to_nbors = generate_neighbor_graph(model)
 
-  n_color = 0
-  cell_color = zeros(Int16, n_cells)
   color_to_inout = Int8[]
+  cell_to_color = zeros(Int16, n_cells)
   touched  = falses(n_cells)
-  q, q_cut = Queue{Int}(), Queue{Int}()
 
-  # First pass: Color IN/OUT cells
-  #   - We assume that every IN/OUT transition can be bridged by a CUT cell
-  first_cell = findfirst(state -> state != CUT, cell_to_inoutcut)
-  enqueue!(q,first_cell); touched[first_cell] = true; # Queue first cell
-  while !isempty(q)
-    cell  = dequeue!(q)
-    nbors = cell_to_nbors[cell]
-    state = cell_to_inoutcut[cell]
-
-    # Mark with color
-    if state != CUT
-      i = findfirst(!iszero,view(cell_color,nbors))
-      if isnothing(i) # New patch
-        n_color += 1
-        cell_color[cell] = n_color
-        push!(color_to_inout, state)
-      else # Existing patch
-        color = cell_color[nbors[i]]
-        cell_color[cell] = color
-      end
-    end
-
-    # Queue and touch unseen neighbors
-    # We touch neighbors here to avoid enqueuing the same cell multiple times
-    for nbor in nbors
-      if !touched[nbor]
-        touched[nbor] = true
-        enqueue!(q,nbor)
-        if cell_to_inoutcut[nbor] == CUT
-          enqueue_interface!(q_cut,cell_to_nbors,cell_to_inoutcut,touched,nbor)
-        end
-      end
+  cell = 1; n_color = zero(Int16)
+  while cell <= n_cells
+    if !touched[cell]
+      n_color += one(Int16)
+      push!(color_to_inout, cell_to_inoutcut[cell])
+      tag_volume!(cell, n_color, cell_to_nbors, cell_to_inoutcut, cell_to_color, touched)
     end
+    cell += 1
   end
 
-  # Second pass: Color CUT cells
-  #   - We assume that every CUT cell has an IN neighbor
-  #   - We assume all IN neighbors have the same color
-  # Then we assign the same color to the CUT cell
-  while !isempty(q_cut)
-    cell  = dequeue!(q_cut)
-    nbors = cell_to_nbors[cell]
-    @assert cell_to_inoutcut[cell] == CUT
+  return cell_to_color, color_to_inout
+end
 
-    i = findfirst(n -> cell_to_inoutcut[n] == IN, nbors)
-    @assert !isnothing(i)
-    cell_color[cell] = cell_color[nbors[i]]
+function find_tagged_volumes(
+  model::DiscreteModel{Dc},tags,
+  cell_to_color::Vector{Int16},
+  color_to_inout::Vector{Int8};
+  Df = Dc - 1
+) where Dc
+  topo = get_grid_topology(model)
+  faces = findall(get_face_mask(get_face_labeling(model),tags,Df))
+  face_to_cell = Geometry.get_faces(topo, Df, Dc)
+
+  is_tagged = falses(length(color_to_inout))
+  for face in faces
+    for cell in view(face_to_cell,face)
+      color = cell_to_color[cell]
+      is_tagged[color] = true
+    end
   end
 
-  return cell_color, color_to_inout
+  return is_tagged
 end
 
 path="./results/FCM_thermal_compliance_ALM_with_islands/"
@@ -206,21 +191,34 @@ for (it,uh,φh,state) in optimiser
     geo = DiscreteGeometry(φh,model)
     bgcell_to_inoutcut = compute_bgcell_to_inoutcut(model,geo)
     colors, color_to_inout = tag_isolated_volumes(model,bgcell_to_inoutcut)
+    color_to_tagged = find_tagged_volumes(model,["Gamma_D"],colors,color_to_inout)
+    cell_to_tagged = map(c -> color_to_tagged[c], colors)
 
     writevtk(Ω,path*"Omega$it",
       cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"velh"=>FEFunction(V_φ,state.vel)],
-      celldata=["inoutcut"=>bgcell_to_inoutcut,"volumes"=>colors];
-      append=false)
+      celldata=[
+        "inoutcut"=>bgcell_to_inoutcut,
+        "volumes"=>colors,
+        "tagged"=>cell_to_tagged,
+        "χ"=>@. (bgcell_to_inoutcut < 0) * (1-cell_to_tagged)
+      ]; append=false)
     writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh])
   end
   write_history(path*"/history.txt",optimiser.history)
+  it == 12 && break
 end
 it = get_history(optimiser).niter; uh = get_state(pcfs)
 geo = DiscreteGeometry(φh,model)
 bgcell_to_inoutcut = compute_bgcell_to_inoutcut(model,geo)
 colors, color_to_inout = tag_isolated_volumes(model,bgcell_to_inoutcut)
+color_to_tagged = find_tagged_volumes(model,["Gamma_D"],colors,color_to_inout)
+cell_to_tagged = map(c -> color_to_tagged[c], colors)
 writevtk(Ω,path*"Omega$it",
   cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh],
-  celldata=["inoutcut"=>bgcell_to_inoutcut,"volumes"=>colors];
-      append=false)
+  celldata=[
+    "inoutcut"=>bgcell_to_inoutcut,
+    "volumes"=>colors,
+    "tagged"=>cell_to_tagged,
+    "χ"=>@. (bgcell_to_inoutcut < 0) * (1-cell_to_tagged)
+  ]; append=false)
 writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh])

scripts/Embedded/Examples/fsi/2-cutfem_navier-stokes_fsi.jl

@@ -2,6 +2,8 @@ using Gridap, Gridap.Geometry, Gridap.Adaptivity
 using GridapEmbedded, GridapEmbedded.LevelSetCutters
 using GridapTopOpt
 
+using GridapSolvers
+
 path = "./results/fsi testing/"
 mkpath(path)
 
@@ -58,7 +60,7 @@ dΓi = Measure(Γi,degree)
 
 # Setup FESpace
 
-uin(x) = VectorValue(0.01x[2]*(1-x[2]),0.0)
+uin(x) = VectorValue(x[2],0.0)
 
 reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P)
 reffe_p = ReferenceFE(lagrangian,Float64,order,space=:P)
@@ -68,7 +70,7 @@ V = TestFESpace(Ω_act,reffe_u,conformity=:H1,dirichlet_tags=["Gamma_D","Gamma_N
 Q = TestFESpace(Ω_act,reffe_p,conformity=:H1)
 T = TestFESpace(Ω_act_solid,reffe_d,conformity=:H1,dirichlet_tags=["Gamma_NoSlipBottom"])
 
-U = TrialFESpace(V,[x->uin(x),VectorValue(0.0,0.0),VectorValue(0.0,0.0)])
+U = TrialFESpace(V,[uin,VectorValue(0.0,0.0),VectorValue(0.0,0.0)])
 P = TrialFESpace(Q)
 R = TrialFESpace(T)
 
@@ -78,24 +80,27 @@ Y = MultiFieldFESpace([V,Q,T])
 # Weak form
 ## Fluid
 # Properties
-μ = 1.0
-ρ = 10.0
+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
 β0 = 0.25
 β1 = 0.2
-β2 = 0.1
-β3 = 0.05
-γ = 10.0
-# Terms
-σf(u,p) = ∇(u)⋅n_Γ - p*n_Γ
+β2 = 0.1*h # 0.05*μ*h
+β3 = 0.05*(μ/h + ρ*u0_max/6)^-1*h^2 # 0.05*h^3
+γ = 100.0
+# Terms (note I've removed 'h's below because they are included in the stabilization parameters)
+σf_n(u,p) = μ*∇(u)⋅n_Γ - p*n_Γ
 a_Ω(u,v) = μ*(∇(u) ⊙ ∇(v))
 b_Ω(v,p) = - (∇⋅v)*p
-c_Γi(p,q) = (β0*h)*jump(p)*jump(q)
+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)
+i_Γg(u,v) = β2*jump(n_Γg⋅∇(u))⋅jump(n_Γg⋅∇(v))
+j_Γg(p,q) = β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Ω +
@@ -112,31 +117,34 @@ jac_fluid((u,p),(du,dp),(v,q)) = a_fluid((du,dp),(v,q)) + dc(u,du,v)
 
 ## Structure
 # Stabilization and material parameters
-γg = 0.1
 function lame_parameters(E,ν)
   λ = (E*ν)/((1+ν)*(1-2*ν))
   μ = E/(2*(1+ν))
   (λ, μ)
 end
-λ, μ = lame_parameters(1.0,0.3)
+λs, μs = lame_parameters(1.0,0.3)
+γg = (λs + 2μs)*0.1
 # Terms
-σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε
+σ(ε) = λs*tr(ε)*one(ε) + 2*μs*ε
 a_solid(d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))dΩout +
-  ∫((γg*h)*jump(n_Γg⋅∇(s)) ⋅ jump(n_Γg⋅∇(d)))dΓg
+  ∫((γg*h^3)*jump(n_Γg⋅∇(s)) ⋅ jump(n_Γg⋅∇(d)))dΓg
 
 ## Full problem
 res((u,p,d),(v,q,s)) = res_fluid((u,p),(v,q)) + a_solid(d,s) +
-  ∫(((σ ∘ ε(d))⋅n_Γ + σf(u,p))⋅s)dΓ # plus sign because of the normal direction
+  ∫(σf_n(u,p) ⋅ s)dΓ # plus sign because of the normal direction
 jac((u,p,d),(du,dp,dd),(v,q,s)) = jac_fluid((u,p),(du,dp),(v,q)) + a_solid(dd,s) +
-  ∫(((σ ∘ ε(dd))⋅n_Γ + σf(du,dp))⋅s)dΓ
+  ∫(σf_n(du,dp) ⋅ s)dΓ
 
 op = FEOperator(res,jac,X,Y)
 solver = GridapSolvers.NewtonSolver(LUSolver();verbose=true)
 uh, ph, dh = solve(solver,op)
 
+# Mass flow rate through surface (this should be close to zero)
+@show m = sum(∫(ρ*uh⋅n_Γ)dΓ)
+
 writevtk(Ω,path*"fsi-navier-stokes-CutFEM_fluid",
   cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh])
 writevtk(Ωout,path*"fsi-navier-stokes-CutFEM_solid",
   cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh])
 
-writevtk(Γ,path*"fsi-navier-stokes-CutFEM_interface",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅n_Γ,"σf"=>σf(uh,ph)])
+writevtk(Γ,path*"fsi-navier-stokes-CutFEM_interface",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅n_Γ,"σf_n"=>σf_n(uh,ph)])

scripts/Embedded/Examples/fsi/2-cutfem_stokes_fsi.jl

@@ -58,7 +58,7 @@ dΓi = Measure(Γi,degree)
 
 # Setup FESpace
 
-uin(x) = VectorValue(0.01x[2]*(1-x[2]),0.0)
+uin(x) = VectorValue(x[2],0.0)
 
 reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P)
 reffe_p = ReferenceFE(lagrangian,Float64,order,space=:P)
@@ -68,7 +68,7 @@ V = TestFESpace(Ω_act,reffe_u,conformity=:H1,dirichlet_tags=["Gamma_D","Gamma_N
 Q = TestFESpace(Ω_act,reffe_p,conformity=:H1)
 T = TestFESpace(Ω_act_solid,reffe_d,conformity=:H1,dirichlet_tags=["Gamma_NoSlipBottom"])
 
-U = TrialFESpace(V,[x->uin(x),VectorValue(0.0,0.0),VectorValue(0.0,0.0)])
+U = TrialFESpace(V,[uin,VectorValue(0.0,0.0),VectorValue(0.0,0.0)])
 P = TrialFESpace(Q)
 R = TrialFESpace(T)
 
@@ -78,23 +78,27 @@ Y = MultiFieldFESpace([V,Q,T])
 # Weak form
 ## Fluid
 # Properties
-μ = 1.0
+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
 β0 = 0.25
 β1 = 0.2
-β2 = 0.1
-β3 = 0.05
-γ = 10.0
+β2 = 0.1*h # 0.05*μ*h
+β3 = 0.05*(μ/h + ρ*u0_max/6)^-1*h^2 # 0.05*h^3
+γ = 100.0
 # Terms
-σf(u,p) = ∇(u)⋅n_Γ - p*n_Γ
+σf_n(u,p) = μ*∇(u)⋅n_Γ - p*n_Γ
 a_Ω(u,v) = μ*(∇(u) ⊙ ∇(v))
 b_Ω(v,p) = - (∇⋅v)*p
-c_Γi(p,q) = (β0*h)*jump(p)*jump(q)
+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)
+i_Γg(u,v) = β2*jump(n_Γg⋅∇(u))⋅jump(n_Γg⋅∇(v))
+j_Γg(p,q) = β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Ω +
@@ -103,30 +107,33 @@ a_fluid((u,p),(v,q)) =
 
 ## Structure
 # Stabilization and material parameters
-γg = 0.1
 function lame_parameters(E,ν)
   λ = (E*ν)/((1+ν)*(1-2*ν))
   μ = E/(2*(1+ν))
   (λ, μ)
 end
-λ, μ = lame_parameters(1.0,0.3)
+λs, μs = lame_parameters(1.0,0.3)
+γg = (λs + 2μs)*0.1
 # Terms
-σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε
+σ(ε) = λs*tr(ε)*one(ε) + 2*μs*ε
 a_solid(d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))dΩout +
-  ∫((γg*h)*jump(n_Γg⋅∇(s)) ⋅ jump(n_Γg⋅∇(d)))dΓg
+  ∫((γg*h^3)*jump(n_Γg⋅∇(s)) ⋅ jump(n_Γg⋅∇(d)))dΓg
 
 ## Full problem
 a((u,p,d),(v,q,s)) = a_fluid((u,p),(v,q)) + a_solid(d,s) +
-  ∫(((σ ∘ ε(d))⋅n_Γ + σf(u,p))⋅s)dΓ # plus sign because of the normal direction
+  ∫(σ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(Ω,path*"fsi-stokes-CutFEM_fluid",
   cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh])
 writevtk(Ωout,path*"fsi-stokes-CutFEM_solid",
   cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh])
 
-writevtk(Γ,path*"fsi-stokes-CutFEM_interface",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅n_Γ,"σf"=>σf(uh,ph)])
+writevtk(Γ,path*"fsi-stokes-CutFEM_interface",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅n_Γ,"σf_n"=>σf_n(uh,ph)])

scripts/Embedded/Examples/stokes & navier-stokes testing/6-Brinkmann_naiver-stokes_P2-P1.jl

@@ -71,7 +71,7 @@ reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P)
 reffe_p = ReferenceFE(lagrangian,Float64,order-1,space=:P)
 
 V = TestFESpace(Ω_act,reffe_u,conformity=:H1,dirichlet_tags=["Gamma_D","Gamma_NoSlip"])
-Q = TestFESpace(Ω_act,reffe_p,conformity=:H1)
+Q = TestFESpace(Ω_act,reffe_p,conformity=:C0)
 
 U = TrialFESpace(V,[x->uin(x),VectorValue(0.0,0.0)])
 P = TrialFESpace(Q)

scripts/Embedded/Examples/stokes & navier-stokes testing/6-Brinkmann_stokes_P2-P1.jl

@@ -66,7 +66,7 @@ reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P)
 reffe_p = ReferenceFE(lagrangian,Float64,order-1,space=:P)
 
 V = TestFESpace(Ω_act,reffe_u,conformity=:H1,dirichlet_tags=["Gamma_D","Gamma_NoSlip"])
-Q = TestFESpace(Ω_act,reffe_p,conformity=:H1)
+Q = TestFESpace(Ω_act,reffe_p,conformity=:C0)
 
 U = TrialFESpace(V,[x->uin(x),VectorValue(0.0,0.0)])
 P = TrialFESpace(Q)