From 892642d188de763cb403407ef85670b658e9a7a5 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Tue, 10 Dec 2024 16:19:44 +1000 Subject: [PATCH] Change to circumdiameter for general h --- docs/src/reference/utilities.md | 2 +- ...Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl | 8 +- ...Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl | 9 +- src/GridapTopOpt.jl | 2 +- .../UnfittedEvolution/CutFEMEvolve.jl | 9 +- .../UnfittedEvolution/StabilisedReinit.jl | 2 +- .../UnfittedEvolution/UnfittedEvolution.jl | 4 +- src/Utilities.jl | 83 ++++++++++++------- test/seq/ElementDiameterTests.jl | 40 +++++++++ 9 files changed, 116 insertions(+), 43 deletions(-) create mode 100644 test/seq/ElementDiameterTests.jl diff --git a/docs/src/reference/utilities.md b/docs/src/reference/utilities.md index 53b1f46..c69a55f 100644 --- a/docs/src/reference/utilities.md +++ b/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 ``` \ No newline at end of file diff --git a/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl b/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl index 5fd79b0..e0112ef 100644 --- a/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-5-Brinkmann_stokes_P1-P1_Ersatz_elast_fsi.jl +++ b/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) diff --git a/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl b/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl index f369798..66a4415 100644 --- a/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl +++ b/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) diff --git a/src/GridapTopOpt.jl b/src/GridapTopOpt.jl index e4899dc..6a298ac 100644 --- a/src/GridapTopOpt.jl +++ b/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") diff --git a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl index 3e4805b..2b637c5 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl +++ b/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: diff --git a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl index 17569e7..f86f9c1 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl +++ b/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) diff --git a/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl b/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl index db8c4af..71b9abf 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl +++ b/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) diff --git a/src/Utilities.jl b/src/Utilities.jl index 4fa14c9..4ce4660 100644 --- a/src/Utilities.jl +++ b/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 \ No newline at end of file diff --git a/test/seq/ElementDiameterTests.jl b/test/seq/ElementDiameterTests.jl new file mode 100644 index 0000000..004c5d4 --- /dev/null +++ b/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 \ No newline at end of file