From 3b483fbac6a8fa7dbfb5dfcfd57409e07c0f7094 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Thu, 3 Oct 2024 17:18:30 +1000 Subject: [PATCH] Bug fix in triangulation and compute derivs --- .../Embedded/SubFacetSkeletonTriangulation.jl | 164 ++++++------------ 1 file changed, 56 insertions(+), 108 deletions(-) diff --git a/scripts/Embedded/SubFacetSkeletonTriangulation.jl b/scripts/Embedded/SubFacetSkeletonTriangulation.jl index ec95f41..27a4e28 100644 --- a/scripts/Embedded/SubFacetSkeletonTriangulation.jl +++ b/scripts/Embedded/SubFacetSkeletonTriangulation.jl @@ -1,5 +1,5 @@ using Gridap -using Gridap.Geometry, Gridap.Adaptivity, Gridap.Algebra, Gridap.Arrays, +using Gridap.Geometry, Gridap.Adaptivity, Gridap.Algebra, Gridap.Arrays, Gridap.FESpaces Gridap.CellData, Gridap.Fields, Gridap.Helpers, Gridap.ReferenceFEs, Gridap.Polynomials using GridapEmbedded @@ -8,6 +8,8 @@ using GridapEmbedded.LevelSetCutters using GridapEmbedded.Interfaces: SubFacetData, SubFacetTriangulation +include("differentiable_trians.jl") + Base.round(a::VectorValue{D,T};kwargs...) where {D,T} = VectorValue{D,T}(round.(a.data;kwargs...)) function SubFacetSkeletonTriangulation(trian::SubFacetTriangulation{Dc,Dp}) where {Dc,Dp} @@ -15,8 +17,9 @@ function SubFacetSkeletonTriangulation(trian::SubFacetTriangulation{Dc,Dp}) wher if Dp == 2 Γ_facet_to_points = map(Reindex(subfacets.point_to_coords),subfacets.facet_to_points) - Γ_pts = unique(round.(reduce(vcat,Γ_facet_to_points);sigdigits=15)) - Γ_facet_to_pts = convert(Vector{Vector{Int32}},map(v->indexin(round.(v;sigdigits=15),Γ_pts),Γ_facet_to_points)) + Γ_pts = unique(x -> round(x;sigdigits=12),reduce(vcat,Γ_facet_to_points)) + Γ_facet_to_pts = convert(Vector{Vector{Int32}},map(v->indexin(round.(v;sigdigits=12), + round.(Γ_pts;sigdigits=12)),Γ_facet_to_points)) grid_top = GridTopology(Γ.subgrid,Table(Γ_facet_to_pts),IdentityVector(length(Γ_pts))) grid = UnstructuredGrid( @@ -61,13 +64,14 @@ end function Γg_to_Γs_perm(Γs::SkeletonTriangulation{0,2},Γg::SkeletonTriangulation{1,2}) grid_top = get_grid_topology(get_background_model(Γg.plus)) + Γ_pts = get_node_coordinates(Γs) # ghost skeleton edges to bg edges Γg_to_bg_edge = Γg.plus.glue.face_to_bgface # bg edges to bg vertices bg_edges_to_bg_vertices = grid_top.n_m_to_nface_to_mfaces[2,1] # points on ghost skeleton Γg_to_bg_vertices = bg_edges_to_bg_vertices[Γg_to_bg_edge] - # points attached to ghost skeleton + # points attached to ghost skeleton Γ_ghost_skel_edges_to_bg_vertex_pts = map(Reindex(grid_top.vertex_coordinates),Γg_to_bg_vertices) # skeleton to vertices Γs_to_vertex = Γs.plus.glue.face_to_bgface @@ -92,23 +96,19 @@ function orient(a::VectorValue{2,T},b::VectorValue{2,T}) where T end order = 1 -n = 10 +n = 51 _model = CartesianDiscreteModel((0,1,0,1),(n,n)) cd = Gridap.Geometry.get_cartesian_descriptor(_model) base_model = UnstructuredDiscreteModel(_model) ref_model = refine(base_model, refinement_method = "barycentric") model = ref_model.model -# model = simplexify(_model) -# model = _model # NOTE: currently broken: -# If we really want this, then need to check that point is on -# boundary of bgcell when computing maps ..._Γ_facet_to_points. Ω = Triangulation(model) dΩ = Measure(Ω,2*order) reffe_scalar = ReferenceFE(lagrangian,Float64,order) V_φ = TestFESpace(model,reffe_scalar) -φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.52,V_φ) # Circle +φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.25,V_φ) # Circle x_φ = get_free_dof_values(φh) if any(isapprox(0.0;atol=10^-10),x_φ) @@ -134,14 +134,12 @@ n2 = CellField(evaluate(get_facet_normal(Γg.minus),Point(0))[Γg_to_Γs],Γs.mi # n1 = GenericCellField(get_facet_normal(Γg.plus),Γs.plus,ReferenceDomain()) # This breaks # n2 = GenericCellField(get_facet_normal(Γg.minus),Γs.minus,ReferenceDomain()) -# n_∂Ω_plus = GenericCellField(get_facet_normal(get_triangulation(Γs.plus)),Γs,ReferenceDomain()) -# n_∂Ω_minus = GenericCellField(get_facet_normal(get_triangulation(Γs.minus)),Γs,ReferenceDomain()) _Γ_trian = Triangulation(get_background_model(Γs)) n_∂Ω_plus = change_domain(get_normal_vector(_Γ_trian).plus,_Γ_trian,ReferenceDomain(),Γs.plus,ReferenceDomain()) n_∂Ω_minus = change_domain(get_normal_vector(_Γ_trian).minus,_Γ_trian,ReferenceDomain(),Γs.minus,ReferenceDomain()) nˢ = Operation(v->(-_2d_cross(v)))(n1) -# nˢ = Operation(orient)(nˢ,n_∂Ω_plus) +# nˢ = Operation(orient)(nˢ,n_∂Ω_plus) # As in other scripts, don't enable this! m1 = Operation(v->(_2d_cross(v)))(n_∂Ω_plus) m2 = Operation(v->(-_2d_cross(v)))(n_∂Ω_minus) @@ -153,8 +151,48 @@ w_Γs = move_to_sub_facet_skeleton(get_fe_basis(V_φ),Γ,Γs) jump_fm = fh_Γs.⁺*m.⁺ + fh_Γs.⁻*m.⁻ # I'm not sure about the basis w_Γs as this is a SkeletonPair... -dJ2 = ∫((nˢ ⋅ (jump_fm*jump(w_Γs)/(abs ∘ (nˢ ⋅ ∇φh_Γs.⁺)))))dΓs -get_array(dJ2) +# dJ2 = ∫(w_Γs*(nˢ ⋅ (jump_fm/(abs ∘ (nˢ ⋅ ∇φh_Γs.⁺)))))dΓs # This is what it should be?? +dJ2_plus = ∫(w_Γs.plus/2*(nˢ ⋅ (jump_fm/(abs ∘ (nˢ ⋅ ∇φh_Γs.⁺)))))dΓs +dJ2_minus = ∫(w_Γs.minus/2*(nˢ ⋅ (jump_fm/(abs ∘ (nˢ ⋅ ∇φh_Γs.⁺)))))dΓs + +dJ2_Γg_plus = map(Reindex(get_array(dJ2_plus)),Γg_to_Γs) # reorder to match GhostSkeleton +dJ2_Γg_minus = map(Reindex(get_array(dJ2_minus)),Γg_to_Γs) # reorder to match GhostSkeleton + +dom_contrib_t2 = DomainContribution() +Gridap.CellData.add_contribution!(dom_contrib_t2,Γg.plus,get_array(dJ2_Γg_plus),-) +Gridap.CellData.add_contribution!(dom_contrib_t2,Γg.minus,get_array(dJ2_Γg_minus),-) +djΓ_exp_vec = assemble_vector(dom_contrib_t2,V_φ) + +## AD +diff_Γ = DifferentiableTriangulation(Γ) +dΓ2 = Measure(diff_Γ,2*order) +jΓ(φ) = ∫(1)dΓ2 +djΓ = gradient(jΓ,φh) +djΓ_contrib = DomainContribution() +Gridap.CellData.add_contribution!(djΓ_contrib,diff_Γ.trian,get_array(djΓ),+) +djΓ_vec_out = assemble_vector(djΓ_contrib,V_φ) + +# Ω_cell_dof = get_cell_dof_ids(V_φ,Ω) +# Γg_cell_dof_plus = get_cell_dof_ids(V_φ,Γg.plus) +# Γg_cell_dof_minus = get_cell_dof_ids(V_φ,Γg.minus) + +# # This is probably barking up the wrong tree @Jordi? +# dv = get_fe_basis(V_φ) +# dv_data = get_data(dv) +# nbasis = length(first(dv_data)) +# contrib1 = [zeros(nbasis) for _ ∈ eachindex(dv_data)] +# contrib2 = [zeros(nbasis) for _ ∈ eachindex(dv_data)] +# for i ∈ eachindex(Γg_cell_dof) +# Ω_cell_dof_idx_plus = findfirst(x -> x == Γg_cell_dof_plus[i],Ω_cell_dof) +# Ω_cell_dof_idx_minus = findfirst(x -> x == Γg_cell_dof_minus[i],Ω_cell_dof) +# contrib1[Ω_cell_dof_idx_plus] = dJ2_Γg_plus[i] +# contrib2[Ω_cell_dof_idx_minus] = dJ2_Γg_minus[i] +# end + +# dom_contrib_t2 = DomainContribution() +# Gridap.CellData.add_contribution!(dom_contrib_t2,Ω,contrib1,-) +# Gridap.CellData.add_contribution!(dom_contrib_t2,Ω,contrib2,-) +# djΓ_exp_vec = assemble_vector(dom_contrib_t2,V_φ) writevtk( Γs, @@ -164,108 +202,18 @@ writevtk( "n_∂Ω_plus"=>n_∂Ω_plus,"n_∂Ω_minus"=>n_∂Ω_minus, "ns"=>nˢ, "∇φh_Γs_plus"=>∇φh_Γs.plus,"∇φh_Γs_minus"=>∇φh_Γs.minus, - "m.minus"=>m.minus,"m.plus"=>m.plus, - "data"=>abs ∘ (∇φh_Γs.plus ⋅ nˢ)] + "m.minus"=>m.minus,"m.plus"=>m.plus] ) writevtk( Ω, "results/Background", - cellfields=["φh"=>∇(φh)] + cellfields=["φh"=>∇(φh),"dj"=>FEFunction(V_φ,djΓ_exp_vec),"dj_AD"=>FEFunction(V_φ,djΓ_vec_out)] ) writevtk( Γ, "results/Boundary" ) writevtk( - Γg.plus.trian, + Γg, "results/GhostSkel" -) - - - -# subfacets.facet_to_points -# Γg.plus.trian.tface_to_mface -bg_edges_to_bg_faces = Ω.model.grid_topology.n_m_to_nface_to_mfaces[2,3] -Γ_ghost_skel_edges_to_bg_edges = bg_edges_to_bg_faces[Γg.plus.glue.face_to_bgface] -subfacets.facet_to_bgcell -Γ_ghost_skel_edges_to_bg_edges -Γ_ghost_skel_edges_to_subfacets = map(v->indexin(v,subfacets.facet_to_bgcell),Γ_ghost_skel_edges_to_bg_edges) - -subfacets_to_Γ_ghost_skel_edges = map(i->findall(v->i ∈ v, Γ_ghost_skel_edges_to_subfacets),eachindex(subfacets.facet_to_bgcell)) -findall(v->isone(length(v)),subfacets_to_Γ_ghost_skel_edges) - -bg_edges_to_bg_faces = Ω.model.grid_topology.n_m_to_nface_to_mfaces[2,1] -Γ_ghost_skel_edges_to_bg_vertices = bg_edges_to_bg_faces[Γg.plus.glue.face_to_bgface] -Γ_ghost_skel_edges_to_bg_vertex_pts = getindex.((Ω.model.grid_topology.vertex_coordinates,),Γ_ghost_skel_edges_to_bg_vertices) - -# map(v->v[2]-v[1],Γ_ghost_skel_edges_to_bg_vertex_pts) - -# Γ_ghost_skel_edges_to_bg_faces[3] -# Γ.subfacets.facet_to_bgcell -# map(Reindex(Γ_facet_to_pts)) - -# getindex.((Γ_pts,),Γ_facet_to_pts) - -subfacets = Γ.subfacets -Γ_facet_to_points = map(Reindex(subfacets.point_to_coords),subfacets.facet_to_points) -# Γ_facet_to_points = sort.(Γ_facet_to_points) - -Γ_pts = unique(round.(reduce(vcat,Γ_facet_to_points);sigdigits=15)) - -Γ_facet_to_pts = convert(Vector{Vector{Int32}},map(v->indexin(round.(v;sigdigits=15),Γ_pts),Γ_facet_to_points)) - -grid_top = GridTopology(Γ.subgrid,Table(Γ_facet_to_pts),IdentityVector(length(Γ_pts))) -grid = UnstructuredGrid( - Γ_pts,Table(Γ_facet_to_pts), - Γ.subgrid.reffes, - Γ.subgrid.cell_types, - Γ.subgrid.orientation_style, - get_facet_normal(Γ) -) - -Γ_model = UnstructuredDiscreteModel(grid,grid_top,FaceLabeling(grid_top)) -Γs = SkeletonTriangulation(Γ_model) - - -# We need a map from points on ghost skeleton to points on skeleton as -# there don't have the same ordering -Γ_pts[Γs.plus.glue.face_to_bgface] # points on skeleton -Γ_ghost_skel_edges_to_bg_vertex_pts # points attached to ghost skeleton - -function belongs_to_edge(q::Point{2,T},p::Vector{Point{2,T}}) where T - p1, p2 = p - tol = 10eps() - @assert abs(p2[1] - p1[1]) > tol || abs(p2[2] - p1[2]) > tol - - if abs(p2[1] - p1[1]) > tol - t = (q[1]-p1[1])/((p2[1]-p1[1])) - 0 < t < 1 && abs(q[2] - (p1[2] + t*(p2[2]-p1[2]))) < tol - else - t = (q[2]-p1[2])/((p2[2]-p1[2])) - 0 < t < 1 && abs(q[1] - (p1[1] + t*(p2[1]-p1[1]))) < tol - end -end - -idxperm = zeros(Int32,length(Γs.plus.glue.face_to_bgface)) -for (i,p) ∈ enumerate(Γ_pts[Γs.plus.glue.face_to_bgface]) - for (j,q) ∈ enumerate(Γ_ghost_skel_edges_to_bg_vertex_pts) - if belongs_to_edge(p,q) - idxperm[i] = j - continue - end - end -end -idxperm - - - -# idxperm = findall.(map(q->belongs_to_edge.(q,Γ_ghost_skel_edges_to_bg_vertex_pts),Γ_pts[Γs.plus.glue.face_to_bgface])) -# @assert all(length.(idxperm) .== 1) -# idxperm = first.(idxperm) - - -writevtk( - Γs, - "results/GammaSkel", - cellfields=["n.plus"=>get_normal_vector(Γs).plus,"n.minus"=>get_normal_vector(Γs).minus] ) \ No newline at end of file