Change to circumdiameter for general h

docs/src/reference/utilities.md

@@ -16,5 +16,5 @@ GridapTopOpt.update_labels!
 GridapTopOpt.initial_lsf
 GridapTopOpt.isotropic_elast_tensor
 GridapTopOpt.get_cartesian_element_sizes
-GridapTopOpt.get_element_sizes
+GridapTopOpt.get_element_diameters
 ```

scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl

@@ -28,8 +28,8 @@ model = GmshDiscreteModel((@__DIR__)*"/mesh.msh")
 writevtk(model,path*"model")
 
 Ω_act = Triangulation(model)
-hₕ = CellField(get_element_sizes(model),Ω_act)
-hmin = GridapTopOpt._get_minimum_element_diameter(hₕ)
+hₕ = CellField(get_element_diameters(model),Ω_act)
+hmin = minimum(get_element_diameters(model))
 
 # Cut the background model
 reffe_scalar = ReferenceFE(lagrangian,Float64,1)
@@ -155,7 +155,7 @@ 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)
+evo = CutFEMEvolve(V_φ,Ω,dΩ_act,hₕ;max_steps,γg=0.1)
 reinit = StabilisedReinit(V_φ,Ω,dΩ_act,hₕ;stabilisation_method=ArtificialViscosity(3.0))
 ls_evo = UnfittedFEEvolution(evo,reinit)
 
@@ -167,7 +167,7 @@ vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
 ## Optimiser
 converged(m) = GridapTopOpt.default_al_converged(
   m;
-  L_tol = 0.5*hmin,
+  L_tol = 0.5hmin,
   C_tol = 0.01vf
 )
 function has_oscillations(m,os_it)

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

@@ -28,7 +28,8 @@ model = GmshDiscreteModel((@__DIR__)*"/mesh.msh")
 writevtk(model,path*"model")
 
 Ω_act = Triangulation(model)
-hₕ = CellField(get_element_sizes(model),Ω_act)
+hₕ = CellField(get_element_diameters(model),Ω_act)
+hmin = minimum(get_element_diameters(model))
 
 # Cut the background model
 reffe_scalar = ReferenceFE(lagrangian,Float64,1)
@@ -107,7 +108,7 @@ cl = H # Characteristic length
 u0_max = maximum(abs,get_dirichlet_dof_values(U))
 μ = ρ*cl*u0_max/Re # Viscosity
 # Stabilization parameters
-γ = 1000.0
+γ(h) = 1e5*μ/h
 
 # Terms
 σf_n(u,p,n) = μ*∇(u) ⋅ n - p*n
@@ -116,7 +117,7 @@ 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
+  ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p) + (γ∘hₕ)*u⋅v ) * Ω.dΩs
 
 ## Structure
 # Stabilization and material parameters
@@ -164,7 +165,7 @@ vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
 ## Optimiser
 converged(m) = GridapTopOpt.default_al_converged(
   m;
-  L_tol = 0.5h,
+  L_tol = 0.5hmin,
   C_tol = 0.01vf
 )
 function has_oscillations(m,os_it)

src/GridapTopOpt.jl

@@ -66,7 +66,7 @@ export update_labels!
 export initial_lsf
 export get_el_Δ
 export get_cartesian_element_sizes
-export get_element_sizes
+export get_element_diameters
 export isotropic_elast_tensor
 
 include("LevelSetEvolution/LevelSetEvolution.jl")

src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl

@@ -47,7 +47,7 @@ get_assembler(s::CutFEMEvolve) = s.assembler
 get_space(s::CutFEMEvolve) = s.space
 get_measure(s::CutFEMEvolve) = s.dΩ_bg
 get_params(s::CutFEMEvolve) = s.params
-get_element_sizes(s::CutFEMEvolve) = s.params.h
+get_element_diameters(s::CutFEMEvolve) = s.params.h
 get_cache(s::CutFEMEvolve) = s.cache
 
 function get_transient_operator(φh,velh,s::CutFEMEvolve)
@@ -58,7 +58,12 @@ function get_transient_operator(φh,velh,s::CutFEMEvolve)
   v_norm = maximum(abs,get_free_dof_values(velh))
   β(vh,∇φ) = vh/(ϵ + v_norm) * ∇φ/(ϵ + norm(∇φ))
 
-  stiffness(t,u,v) = ∫(((β ∘ (velh,∇(φh))) ⋅ ∇(u)) * v)dΩ_bg + ∫((γg*h*h)*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg
+  γ(h) = γg*h^2
+
+  aₛ(u,v,h::CellField) = ∫(mean(γ ∘ h)*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg
+  aₛ(u,v,h::Real) = ∫(γ(h)*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg
+
+  stiffness(t,u,v) = ∫(((β ∘ (velh,∇(φh))) ⋅ ∇(u)) * v)dΩ_bg + aₛ(u,v,h)
   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:

src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl

@@ -64,7 +64,7 @@ get_space(s::StabilisedReinit) = s.space
 get_embedded_collection(s::StabilisedReinit) = s.Ωs
 get_measure(s::StabilisedReinit) = s.dΩ_bg
 get_params(s::StabilisedReinit) = s.params
-get_element_sizes(s::StabilisedReinit) = s.params.h
+get_element_diameters(s::StabilisedReinit) = s.params.h
 get_cache(s::StabilisedReinit) = s.cache
 
 function solve!(s::StabilisedReinit,φh,nls_cache::Nothing)

src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl

@@ -42,8 +42,8 @@ struct UnfittedFEEvolution{A<:Evolver,B<:Reinitialiser,C<:Real} <: LevelSetEvolu
   reinitialiser::B
   min_element_diameter::C
   function UnfittedFEEvolution(evolver::A,reinitialiser::B) where {A,B}
-    h_evo = get_element_sizes(evolver)
-    h_reinit = get_element_sizes(reinitialiser)
+    h_evo = get_element_diameters(evolver)
+    h_reinit = get_element_diameters(reinitialiser)
     @assert h_evo === h_reinit "Element sizes for evolution and reinitialisation should be the same."
     min_element_diameter = _get_minimum_element_diameter(h_evo)
     new{A,B,typeof(min_element_diameter)}(evolver,reinitialiser,min_element_diameter)

src/Utilities.jl

@@ -244,44 +244,71 @@ function get_order(space::DistributedFESpace)
   return getany(order)
 end
 
-# Element size for general model
+# Element diameter for general model
 
+# According to Massing et al. (doi:10.1007/s10915-014-9838-9) the stabilisation terms
+#   should use h_K denotes the diameter of element K, and h_F denotes the average
+#   of the diameters of the elements sharing a facet K. The latter can be computed
+#   as the mean on the facet triangulation.
+#
+# In this implementation, diameter is interpreted as the circumdiameter of the polytope.
 """
-    get_element_sizes(model)
+    get_element_diameters(model)
 
-Given a general model return the element size as a lazy map of minimum (by default)
-side length sizes.
+Given a general unstructured model return the element circumdiameters.
 """
-function get_element_sizes(model,f=min)
+function get_element_diameters(model)
   coords = get_cell_coordinates(model)
   polys = get_polytopes(model)
   @assert length(polys) == 1 "Only one cell type is currently supported"
-  return _get_element_sizes(f,coords,first(polys))
-end
-
-function _get_tri_side_lengths(cell_coords)
-  p1 = cell_coords[1]
-  p2 = cell_coords[2]
-  p3 = cell_coords[3]
-  return norm(p1-p2), norm(p1-p3), norm(p2-p3)
-end
-
-function _get_tet_side_lengths(cell_coords)
-  p1 = cell_coords[1]
-  p2 = cell_coords[2]
-  p3 = cell_coords[3]
-  p4 = cell_coords[4]
-  return norm(p1-p2), norm(p1-p3), norm(p2-p3),
-    norm(p1-p4), norm(p2-p4), norm(p3-p4)
-end
-
-function _get_element_sizes(f,coords,poly::Polytope)
+  poly = first(polys)
   if poly == TRI
-    side_lengths = lazy_map(_get_tri_side_lengths,coords)
+    return lazy_map(_get_tri_circumdiameter,coords)
   elseif poly == TET
-    side_lengths = lazy_map(_get_tet_side_lengths,coords)
+    return lazy_map(_get_tet_circumdiameter,coords)
   else
     @notimplemented "Only triangles and tetrahedra are currently supported"
   end
-  return lazy_map(t->f(t...),side_lengths)
+end
+
+# Based on doi:10.1017/CBO9780511973611. C is the Cayley-Menger matrix.
+function _get_tri_circumdiameter(coords)
+  d12 = norm(coords[1]-coords[2])^2
+  d13 = norm(coords[1]-coords[3])^2
+  d23 = norm(coords[2]-coords[3])^2
+  # C = [
+  #   0  1   1   1
+  #   1  0  d12 d13
+  #   1 d12  0  d23
+  #   1 d13 d23  0
+  # ];
+  # M = -2inv(C);
+  # circumcentre = (M[1,2]*coords[1] + M[1,3]*coords[2] + M[1,4]*coords[3])/sum(M[1,2:end])
+  # circumdiameter = sqrt(M[1,1])
+  M11 = -((4*d12*d13*d23)/(d12^2+(d13-d23)^2-2*d12*(d13+d23)))
+  return sqrt(M11)
+end
+
+function _get_tet_circumdiameter(coords)
+  d12 = norm(coords[1]-coords[2])^2
+  d13 = norm(coords[1]-coords[3])^2
+  d14 = norm(coords[1]-coords[4])^2
+  d23 = norm(coords[2]-coords[3])^2
+  d24 = norm(coords[2]-coords[4])^2
+  d34 = norm(coords[3]-coords[4])^2
+  # C = [
+  #   0  1   1   1   1
+  #   1  0  d12 d13 d14
+  #   1 d12  0  d23 d24
+  #   1 d13 d23  0  d34
+  #   1 d14 d24 d34  0
+  # ];
+  # M = -2inv(C);
+  # circumcentre = (M[1,2]*coords[1] + M[1,3]*coords[2] + M[1,4]*coords[3] + M[1,5]*coords[4])/sum(M[1,2:end])
+  # circumdiameter = sqrt(M[1,1])
+  M11 = (d14^2*d23^2+(d13*d24-d12*d34)^2-2*d14*d23*(d13*d24+d12*d34))/(
+    d13^2*d24+d12^2*d34+d23*(d14^2+d14*(d23-d24-d34)+d24*d34)-
+    d13*(d14*(d23+d24-d34)+d24*(d23-d24+d34))-d12*((d23+d24-d34)*d34+
+    d14*(d23-d24+d34)+d13*(-d23+d24+d34)))
+  return sqrt(M11)
 end

test/seq/ElementDiameterTests.jl

@@ -0,0 +1,40 @@
+module ElementDiameterTests
+using Test
+using GridapTopOpt
+using GridapTopOpt: _get_tri_circumdiameter, _get_tet_circumdiameter
+
+function generate_model(D,n)
+  domain = (D==2) ? (0,1,0,1) : (0,1,0,1,0,1)
+  cell_partition = (D==2) ? (n,n) : (n,n,n)
+  base_model = UnstructuredDiscreteModel((CartesianDiscreteModel(domain,cell_partition)))
+  ref_model = refine(base_model, refinement_method = "barycentric")
+  model = ref_model.model
+  return model
+end
+
+tri_coords = [
+  VectorValue(0.2948854413908402,0.2331928283531077),
+  VectorValue(0.2886157430114769,0.2306953723644822),
+  VectorValue(0.2908695476768671,0.2387503094465989)
+];
+
+@test abs(_get_tri_circumdiameter(tri_coords)/2 - 0.00431269)/0.0043126 < 1e-5
+
+tet_coords = [
+  VectorValue(0.2,0.2,0.1),
+  VectorValue(1,0.2,0.1),
+  VectorValue(2.3,1,0.1),
+  VectorValue(1.2,3/2,1)
+];
+
+@test abs(_get_tet_circumdiameter(tet_coords)/2 - 2.64105)/2.64105 < 1e-5
+
+model = generate_model(2,1)
+@test all(isone,get_element_diameters(model))
+
+get_cell_coordinates(model)
+
+model = generate_model(3,1)
+@test all(abs.(get_element_diameters(model)/2 .- 0.559017) .< 1e-8)
+
+end