From 72f10d07367a3f49550e617969d193b8d40e40cb Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 27 Jan 2023 09:54:38 +0100 Subject: [PATCH 1/9] Minor change --- src/Adaptivity/Adaptivity.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Adaptivity/Adaptivity.jl b/src/Adaptivity/Adaptivity.jl index 128bc32a4..9cc5971a4 100644 --- a/src/Adaptivity/Adaptivity.jl +++ b/src/Adaptivity/Adaptivity.jl @@ -47,7 +47,6 @@ include("AdaptivityGlues.jl") include("AdaptedDiscreteModels.jl") include("AdaptedTriangulations.jl") include("CompositeQuadratures.jl") - include("EdgeBasedRefinement.jl") end # module From d2eb3080d819b3e079a2ac005a7e202064e34ea8 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 9 Feb 2023 11:38:53 +1100 Subject: [PATCH 2/9] adapted change_domain now returns similar_cell_field --- src/Adaptivity/AdaptedTriangulations.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Adaptivity/AdaptedTriangulations.jl b/src/Adaptivity/AdaptedTriangulations.jl index 73337d78e..e75b5dc00 100644 --- a/src/Adaptivity/AdaptedTriangulations.jl +++ b/src/Adaptivity/AdaptedTriangulations.jl @@ -242,10 +242,10 @@ function change_domain_o2n(f_coarse,ftrian::AdaptedTriangulation{Dc},glue::Adapt # Final map: f_fine(x_fine) = f_f2c ∘ Φ^(-1)(x_fine) = f_coarse(x_coarse) f_fine = lazy_map(∘,f_f2c,ref_coord_map) - return GenericCellField(f_fine,ftrian,ReferenceDomain()) + return CellData.similar_cell_field(f_coarse,f_fine,ftrian,ReferenceDomain()) else f_fine = Fill(Gridap.Fields.ConstantField(0.0),num_cells(ftrian)) - return GenericCellField(f_fine,ftrian,ReferenceDomain()) + return CellData.similar_cell_field(f_coarse,f_fine,ftrian,ReferenceDomain()) end end @@ -269,9 +269,9 @@ function change_domain_n2o(f_fine,ctrian::Triangulation{Dc},glue::AdaptivityGlue rrules = get_old_cell_refinement_rules(glue) f_coarse = lazy_map(FineToCoarseField,f_c2f,rrules) - return GenericCellField(f_coarse,ctrian,ReferenceDomain()) + return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain()) else f_coarse = Fill(Gridap.Fields.ConstantField(0.0),num_cells(ftrian)) - return GenericCellField(f_coarse,ctrian,ReferenceDomain()) + return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain()) end end From 288114d22c628225afc0c854905685dedcd10bbc Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 13 Feb 2023 19:33:20 +1100 Subject: [PATCH 3/9] Added adaptivity change_domain support for views and boundary trians --- src/Adaptivity/AdaptedTriangulations.jl | 64 +++++++++---- .../ExtendingChangeDomainTests.jl | 91 +++++++++++++++++++ 2 files changed, 138 insertions(+), 17 deletions(-) create mode 100644 test/AdaptivityTests/ExtendingChangeDomainTests.jl diff --git a/src/Adaptivity/AdaptedTriangulations.jl b/src/Adaptivity/AdaptedTriangulations.jl index e75b5dc00..255c5d590 100644 --- a/src/Adaptivity/AdaptedTriangulations.jl +++ b/src/Adaptivity/AdaptedTriangulations.jl @@ -112,9 +112,9 @@ for (stype,ttype) in [(:AdaptedTriangulation,:AdaptedTriangulation),(:AdaptedTri return is_change_possible($sstrian,$tttrian) end - if isa($sstrian,BodyFittedTriangulation) && isa($tttrian,BodyFittedTriangulation) + #if isa($sstrian,BodyFittedTriangulation) && isa($tttrian,BodyFittedTriangulation) return is_related(strian,ttrian) - end + #end # Support for Boundary, Skeletton, ... triangulations is not yet implemented # unless they are based on the same background model. @@ -151,7 +151,8 @@ function CellData.change_domain(a::CellField,strian::Triangulation,::ReferenceDo @check is_change_possible(strian,ttrian) if (get_background_model(strian) === get_background_model(ttrian)) - return change_domain(a,strian,ReferenceDomain(),ttrian.trian,ReferenceDomain()) + b = change_domain(a,strian,ReferenceDomain(),ttrian.trian,ReferenceDomain()) + return CellData.similar_cell_field(b,get_data(b),ttrian,ReferenceDomain()) end return change_domain_o2n(a,ttrian) @@ -206,15 +207,15 @@ function Geometry.move_contributions(scell_to_val::AbstractArray, glue::Adaptivi return tcell_to_val end -function change_domain_o2n(f_old, new_trian::AdaptedTriangulation{Dc,Dp}) where {Dc,Dp} - @notimplementedif num_dims(get_triangulation(f_old)) != Dc +function change_domain_o2n(f_old,new_trian::AdaptedTriangulation{Dc,Dp}) where {Dc,Dp} + #@notimplementedif num_dims(get_triangulation(f_old)) != Dc glue = get_adaptivity_glue(new_trian) return change_domain_o2n(f_old,new_trian,glue) end -function change_domain_n2o(f_new, old_trian::Triangulation{Dc,Dp}) where {Dc,Dp} +function change_domain_n2o(f_new,old_trian::Triangulation{Dc,Dp}) where {Dc,Dp} new_trian = get_triangulation(f_new) - @notimplementedif num_dims(new_trian) != Dc + #@notimplementedif num_dims(new_trian) != Dc @check isa(new_trian,AdaptedTriangulation) glue = get_adaptivity_glue(new_trian) return change_domain_n2o(f_new,old_trian,glue) @@ -228,20 +229,34 @@ function change_domain_o2n(f_old,new_trian::AdaptedTriangulation,glue::Adaptivit @notimplemented end -function change_domain_o2n(f_coarse,ftrian::AdaptedTriangulation{Dc},glue::AdaptivityGlue{<:RefinementGlue,Dc}) where Dc +function change_domain_o2n(f_coarse,ftrian::AdaptedTriangulation{Dc},glue::AdaptivityGlue{<:RefinementGlue}) where Dc ctrian = get_triangulation(f_coarse) - @notimplementedif num_dims(ctrian) != Dc + #@notimplementedif num_dims(ctrian) != Dc + + D = num_cell_dims(ctrian) + cglue = get_glue(ctrian,Val(D)) + fglue = get_glue(ftrian,Val(D)) if (num_cells(ctrian) != 0) + ### Old Triangulation -> Old Model + coarse_tface_to_field = CellData.get_data(f_coarse) + coarse_mface_to_field = extend(coarse_tface_to_field,cglue.mface_to_tface) + + ### Old Model -> New Model # Coarse field but with fine indexing, i.e # f_f2c[i_fine] = f_coarse[coarse_parent(i_fine)] - f_f2c = c2f_reindex(f_coarse,glue) + f_f2c = c2f_reindex(coarse_mface_to_field,glue) # Fine to coarse coordinate map: x_coarse = Φ^(-1)(x_fine) ref_coord_map = get_n2o_reference_coordinate_map(glue) - # Final map: f_fine(x_fine) = f_f2c ∘ Φ^(-1)(x_fine) = f_coarse(x_coarse) - f_fine = lazy_map(∘,f_f2c,ref_coord_map) + # f_fine(x_fine) = f_f2c ∘ Φ^(-1)(x_fine) = f_coarse(x_coarse) + fine_mface_to_field = lazy_map(∘,f_f2c,ref_coord_map) + + ### New Model -> New Triangulation + fine_tface_to_field = lazy_map(Reindex(fine_mface_to_field),fglue.tface_to_mface) + f_fine = lazy_map(Broadcasting(∘),fine_tface_to_field,fglue.tface_to_mface_map) + return CellData.similar_cell_field(f_coarse,f_fine,ftrian,ReferenceDomain()) else f_fine = Fill(Gridap.Fields.ConstantField(0.0),num_cells(ftrian)) @@ -257,21 +272,36 @@ function change_domain_n2o(f_new,old_trian::Triangulation,glue::AdaptivityGlue) @notimplemented end -function change_domain_n2o(f_fine,ctrian::Triangulation{Dc},glue::AdaptivityGlue{<:RefinementGlue,Dc}) where Dc - @notimplementedif num_dims(ctrian) != Dc +function change_domain_n2o(f_fine,ctrian::Triangulation{Dc},glue::AdaptivityGlue{<:RefinementGlue}) where Dc + ftrian = get_triangulation(f_fine) + #@notimplementedif num_dims(ftrian) != Dc msg = "Evaluating a fine CellField in the coarse mesh is costly! If you are using this feature to integrate, consider using a CompositeMeasure instead (see test/AdaptivityTests/GridTransferTests.jl)." @warn msg + D = num_cell_dims(ftrian) + cglue = get_glue(ctrian,Val(D)) + fglue = get_glue(ftrian,Val(D)) + if (num_cells(ctrian) != 0) + ### New Triangulation -> New Model + fine_tface_to_field = CellData.get_data(f_fine) + fine_mface_to_field = extend(fine_tface_to_field,fglue.mface_to_tface) + + ### New Model -> Old Model # f_c2f[i_coarse] = [f_fine[i_fine_1], ..., f_fine[i_fine_nChildren]] - f_c2f = f2c_reindex(f_fine,glue) + f_c2f = f2c_reindex(fine_mface_to_field,glue) rrules = get_old_cell_refinement_rules(glue) - f_coarse = lazy_map(FineToCoarseField,f_c2f,rrules) + coarse_mface_to_field = lazy_map(FineToCoarseField,f_c2f,rrules) + + ### Old Model -> Old Triangulation + coarse_tface_to_field = lazy_map(Reindex(coarse_mface_to_field),cglue.tface_to_mface) + f_coarse = lazy_map(Broadcasting(∘),coarse_tface_to_field,cglue.tface_to_mface_map) + return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain()) else - f_coarse = Fill(Gridap.Fields.ConstantField(0.0),num_cells(ftrian)) + f_coarse = Fill(Gridap.Fields.ConstantField(0.0),num_cells(fcoarse)) return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain()) end end diff --git a/test/AdaptivityTests/ExtendingChangeDomainTests.jl b/test/AdaptivityTests/ExtendingChangeDomainTests.jl new file mode 100644 index 000000000..b4e0828db --- /dev/null +++ b/test/AdaptivityTests/ExtendingChangeDomainTests.jl @@ -0,0 +1,91 @@ +module ExtendingChangeDomainTests + +using Test +using Gridap +using Gridap.Geometry +using Gridap.CellData +using Gridap.Adaptivity +using Gridap.ReferenceFEs +using FillArrays + +norm2(u,Ω::Triangulation) = sum(∫(u⋅u)*Measure(Ω,qorder)) +norm2(u,dΩ::Measure) = sum(∫(u⋅u)*dΩ) + +order = 1 +qorder = 2*order+1 +sol(x) = sum(x) + +D = 2 +domain = Tuple(repeat([0,1],D)) +cmodel = CartesianDiscreteModel(domain,Tuple(fill(2,D))) +fmodel = refine(cmodel,2) + +Ωc = Triangulation(cmodel) +Ωf = Triangulation(fmodel) + +reffe = ReferenceFE(lagrangian,Float64,order) +Vc = TestFESpace(cmodel,reffe) +Uc = TrialFESpace(Vc) +Vf = TestFESpace(fmodel,reffe) +Uf = TrialFESpace(Vf) + +""" + BodyFittedTriangulation Views +""" + +Ωc_view = view(Ωc,[1,2]) +Ωf_view = view(Ωf,[1,2,3,4,5,6,7,8]) + +v_Ωf = interpolate(sol,Vf) +v_Ωf_view = change_domain(v_Ωf,Ωf_view,ReferenceDomain()) + +v_Ωc_view_view = change_domain(v_Ωf_view,Ωc_view,ReferenceDomain()) + +v_Ωc = interpolate(sol,Vc) +v_Ωc_view = change_domain(v_Ωc,Ωc_view,ReferenceDomain()) + +v_Ωf_view_view = change_domain(v_Ωc_view,Ωf_view,ReferenceDomain()) + +@test norm2(v_Ωc_view_view,Ωc_view) ≈ norm2(v_Ωc_view,Ωc_view) +@test norm2(v_Ωf_view_view,Ωf_view) ≈ norm2(v_Ωf_view,Ωf_view) + +# Same but automatic +@test norm2(v_Ωc_view,Ωf_view) ≈ norm2(v_Ωf_view,Ωf_view) +@test norm2(v_Ωf_view,Ωc_view) ≈ norm2(v_Ωc_view,Ωc_view) + +""" + BodyFitted -> Boundary +""" +Γc = Boundary(cmodel) +Γf = Boundary(fmodel) + +v_Ωc = interpolate(sol,Vc) +v_Ωf = interpolate(sol,Vf) + +v_Γf = change_domain(v_Ωf,Γf,ReferenceDomain()) +v_Γc = change_domain(v_Ωc,Γc,ReferenceDomain()) + +@test norm2(v_Ωf,Γc) ≈ norm2(v_Γc,Γc) +@test norm2(v_Ωc,Γf) ≈ norm2(v_Γf,Γf) + +""" + BodyFitted -> Skeleton + + Not currently working, we need to treat SkeletonPairs separately.... +""" + +""" +Λc = Skeleton(cmodel) +Λf = Skeleton(fmodel) + +v_Ωc = interpolate(sol,Vc).plus +v_Ωf = interpolate(sol,Vf).plus + +v_Λf = change_domain(v_Ωf,Λf,ReferenceDomain()) +v_Λc = change_domain(v_Ωc,Λc,ReferenceDomain()) + +@test norm2(v_Ωf,Λc) ≈ norm2(v_Λc,Λc) +@test norm2(v_Ωc,Λf) ≈ norm2(v_Λf,Λf) +""" + +end \ No newline at end of file From 5ca897dd5dd116856220058ae99c5e2dfc5cb4df Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 14 Feb 2023 09:37:35 +1100 Subject: [PATCH 4/9] Added adaptivity change_domain support for SkeletonTriangulations --- src/Adaptivity/AdaptedTriangulations.jl | 83 ++++++++++++------- .../ExtendingChangeDomainTests.jl | 16 ++-- 2 files changed, 62 insertions(+), 37 deletions(-) diff --git a/src/Adaptivity/AdaptedTriangulations.jl b/src/Adaptivity/AdaptedTriangulations.jl index 255c5d590..a3de3d982 100644 --- a/src/Adaptivity/AdaptedTriangulations.jl +++ b/src/Adaptivity/AdaptedTriangulations.jl @@ -155,7 +155,7 @@ function CellData.change_domain(a::CellField,strian::Triangulation,::ReferenceDo return CellData.similar_cell_field(b,get_data(b),ttrian,ReferenceDomain()) end - return change_domain_o2n(a,ttrian) + return change_domain_o2n(a,strian,ttrian) end function CellData.change_domain(a::CellField,strian::AdaptedTriangulation,::ReferenceDomain,ttrian::Triangulation,::ReferenceDomain) @@ -168,7 +168,7 @@ function CellData.change_domain(a::CellField,strian::AdaptedTriangulation,::Refe return change_domain(a,strian.trian,ReferenceDomain(),ttrian,ReferenceDomain()) end - return change_domain_n2o(a,ttrian) + return change_domain_n2o(a,strian,ttrian) end function CellData.change_domain(a::CellField,strian::AdaptedTriangulation,::ReferenceDomain,ttrian::AdaptedTriangulation,::ReferenceDomain) @@ -207,35 +207,29 @@ function Geometry.move_contributions(scell_to_val::AbstractArray, glue::Adaptivi return tcell_to_val end -function change_domain_o2n(f_old,new_trian::AdaptedTriangulation{Dc,Dp}) where {Dc,Dp} - #@notimplementedif num_dims(get_triangulation(f_old)) != Dc +# Change domain o2n/n2o + +function change_domain_o2n(f_old,old_trian::Triangulation,new_trian::AdaptedTriangulation) glue = get_adaptivity_glue(new_trian) - return change_domain_o2n(f_old,new_trian,glue) + return change_domain_o2n(f_old,old_trian,new_trian,glue) end -function change_domain_n2o(f_new,old_trian::Triangulation{Dc,Dp}) where {Dc,Dp} - new_trian = get_triangulation(f_new) - #@notimplementedif num_dims(new_trian) != Dc - @check isa(new_trian,AdaptedTriangulation) +function change_domain_n2o(f_new,new_trian::AdaptedTriangulation,old_trian::Triangulation) glue = get_adaptivity_glue(new_trian) - return change_domain_n2o(f_new,old_trian,glue) + return change_domain_n2o(f_new,new_trian,old_trian,glue) end """ Given a AdaptivityGlue and a CellField defined on the parent(old) mesh, returns an equivalent CellField on the child(new) mesh. """ -function change_domain_o2n(f_old,new_trian::AdaptedTriangulation,glue::AdaptivityGlue) +function change_domain_o2n(f_old,old_trian::Triangulation,new_trian::AdaptedTriangulation,glue::AdaptivityGlue) @notimplemented end -function change_domain_o2n(f_coarse,ftrian::AdaptedTriangulation{Dc},glue::AdaptivityGlue{<:RefinementGlue}) where Dc - ctrian = get_triangulation(f_coarse) - #@notimplementedif num_dims(ctrian) != Dc - - D = num_cell_dims(ctrian) - cglue = get_glue(ctrian,Val(D)) - fglue = get_glue(ftrian,Val(D)) +function change_domain_o2n(f_coarse,ctrian::Triangulation{Dc},ftrian::AdaptedTriangulation,glue::AdaptivityGlue{<:RefinementGlue}) where Dc + cglue = get_glue(ctrian,Val(Dc)) + fglue = get_glue(ftrian,Val(Dc)) if (num_cells(ctrian) != 0) ### Old Triangulation -> Old Model @@ -268,20 +262,13 @@ end Given a AdaptivityGlue and a CellField defined on the child(new) mesh, returns an equivalent CellField on the parent(old) mesh. """ -function change_domain_n2o(f_new,old_trian::Triangulation,glue::AdaptivityGlue) +function change_domain_n2o(f_new,new_trian::AdaptedTriangulation,old_trian::Triangulation,glue::AdaptivityGlue) @notimplemented end -function change_domain_n2o(f_fine,ctrian::Triangulation{Dc},glue::AdaptivityGlue{<:RefinementGlue}) where Dc - ftrian = get_triangulation(f_fine) - #@notimplementedif num_dims(ftrian) != Dc - msg = "Evaluating a fine CellField in the coarse mesh is costly! If you are using this feature - to integrate, consider using a CompositeMeasure instead (see test/AdaptivityTests/GridTransferTests.jl)." - @warn msg - - D = num_cell_dims(ftrian) - cglue = get_glue(ctrian,Val(D)) - fglue = get_glue(ftrian,Val(D)) +function change_domain_n2o(f_fine,ftrian::AdaptedTriangulation{Dc},ctrian::Triangulation,glue::AdaptivityGlue{<:RefinementGlue}) where Dc + cglue = get_glue(ctrian,Val(Dc)) + fglue = get_glue(ftrian,Val(Dc)) if (num_cells(ctrian) != 0) ### New Triangulation -> New Model @@ -305,3 +292,41 @@ function change_domain_n2o(f_fine,ctrian::Triangulation{Dc},glue::AdaptivityGlue return CellData.similar_cell_field(f_fine,f_coarse,ctrian,ReferenceDomain()) end end + + +# Specialisation for Skeleton Pairs + +function change_domain_o2n(f_old,old_trian::Triangulation,new_trian::AdaptedTriangulation{Dc,Dp,<:SkeletonTriangulation}) where {Dc,Dp} + @check isa(CellData.get_data(f_old),Fill) + new_trian_plus = AdaptedTriangulation(new_trian.trian.plus,new_trian.adapted_model) + return change_domain_o2n(f_old,old_trian,new_trian_plus) +end + +function change_domain_n2o(f_new,new_trian::AdaptedTriangulation,old_trian::SkeletonTriangulation) + @check isa(CellData.get_data(f_new),Fill) + return change_domain_n2o(f_new,new_trian,old_trian.plus) +end + +function change_domain_o2n(f_old::CellData.CellFieldAt,old_trian::Triangulation,new_trian::AdaptedTriangulation{Dc,Dp,<:SkeletonTriangulation}) where {Dc,Dp} + if isa(f_old,CellData.CellFieldAt{:plus}) + _new_trian = AdaptedTriangulation(new_trian.trian.plus,new_trian.adapted_model) + elseif isa(f_old,CellData.CellFieldAt{:minus}) + _new_trian = AdaptedTriangulation(new_trian.trian.minus,new_trian.adapted_model) + else + @unreachable + end + return change_domain_o2n(f_old,old_trian,_new_trian) +end + +function change_domain_n2o(f_new::CellData.CellFieldAt,new_trian::AdaptedTriangulation,old_trian::SkeletonTriangulation) + if isa(f_new,CellData.CellFieldAt{:plus}) + _old_trian = old_trian.plus + elseif isa(f_new,CellData.CellFieldAt{:minus}) + _old_trian = old_trian.minus + else + @unreachable + end + return change_domain_n2o(f_new,new_trian,_old_trian) +end + + diff --git a/test/AdaptivityTests/ExtendingChangeDomainTests.jl b/test/AdaptivityTests/ExtendingChangeDomainTests.jl index b4e0828db..242050a26 100644 --- a/test/AdaptivityTests/ExtendingChangeDomainTests.jl +++ b/test/AdaptivityTests/ExtendingChangeDomainTests.jl @@ -74,18 +74,18 @@ v_Γc = change_domain(v_Ωc,Γc,ReferenceDomain()) Not currently working, we need to treat SkeletonPairs separately.... """ -""" Λc = Skeleton(cmodel) Λf = Skeleton(fmodel) -v_Ωc = interpolate(sol,Vc).plus -v_Ωf = interpolate(sol,Vf).plus +for sign in [:plus,:minus] + v_Ωc = getproperty(interpolate(sol,Vc),sign) + v_Ωf = getproperty(interpolate(sol,Vf),sign) -v_Λf = change_domain(v_Ωf,Λf,ReferenceDomain()) -v_Λc = change_domain(v_Ωc,Λc,ReferenceDomain()) + v_Λf = change_domain(v_Ωf,Λf,ReferenceDomain()) + v_Λc = change_domain(v_Ωc,Λc,ReferenceDomain()) -@test norm2(v_Ωf,Λc) ≈ norm2(v_Λc,Λc) -@test norm2(v_Ωc,Λf) ≈ norm2(v_Λf,Λf) -""" + @test norm2(v_Ωf,Λc) ≈ norm2(v_Λc,Λc) + @test norm2(v_Ωc,Λf) ≈ norm2(v_Λf,Λf) +end end \ No newline at end of file From 8c64ca1e04a7ec91790e8dfe8a513b1420060a8d Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 14 Feb 2023 10:09:43 +1100 Subject: [PATCH 5/9] Added new @notimplemented flags --- src/Adaptivity/AdaptedTriangulations.jl | 35 ++++++++++++------------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/src/Adaptivity/AdaptedTriangulations.jl b/src/Adaptivity/AdaptedTriangulations.jl index a3de3d982..f25e16510 100644 --- a/src/Adaptivity/AdaptedTriangulations.jl +++ b/src/Adaptivity/AdaptedTriangulations.jl @@ -112,14 +112,9 @@ for (stype,ttype) in [(:AdaptedTriangulation,:AdaptedTriangulation),(:AdaptedTri return is_change_possible($sstrian,$tttrian) end - #if isa($sstrian,BodyFittedTriangulation) && isa($tttrian,BodyFittedTriangulation) - return is_related(strian,ttrian) - #end - - # Support for Boundary, Skeletton, ... triangulations is not yet implemented - # unless they are based on the same background model. - @notimplemented - return false + strian_is_cell_wise = (num_cell_dims(strian) == num_point_dims(strian)) + trians_are_related = is_related(strian,ttrian) + return strian_is_cell_wise && trians_are_related end end @eval begin @@ -131,13 +126,17 @@ for (stype,ttype) in [(:AdaptedTriangulation,:AdaptedTriangulation),(:AdaptedTri if (get_background_model(strian) === get_background_model(ttrian)) return best_target($sstrian,$tttrian) end + + ttrian_is_cell_wise = (num_cell_dims(ttrian) == num_point_dims(ttrian)) + strian_is_cell_wise = (num_cell_dims(strian) == num_point_dims(strian)) - if isa($sstrian,BodyFittedTriangulation) && isa($tttrian,BodyFittedTriangulation) + if strian_is_cell_wise && ttrian_is_cell_wise + # Choose the child mesh, both ways are possible is_child(ttrian,strian) ? (return ttrian) : (return strian) + elseif strian_is_cell_wise + return strian end - - # Support for Boundary, Skeletton, ... triangulations is not yet implemented - # unless they are based on the same background model. + @notimplemented return nothing end @@ -145,9 +144,6 @@ for (stype,ttype) in [(:AdaptedTriangulation,:AdaptedTriangulation),(:AdaptedTri end function CellData.change_domain(a::CellField,strian::Triangulation,::ReferenceDomain,ttrian::AdaptedTriangulation,::ReferenceDomain) - if strian === ttrian - return a - end @check is_change_possible(strian,ttrian) if (get_background_model(strian) === get_background_model(ttrian)) @@ -159,9 +155,6 @@ function CellData.change_domain(a::CellField,strian::Triangulation,::ReferenceDo end function CellData.change_domain(a::CellField,strian::AdaptedTriangulation,::ReferenceDomain,ttrian::Triangulation,::ReferenceDomain) - if strian === ttrian - return a - end @check is_change_possible(strian,ttrian) if (get_background_model(strian) === get_background_model(ttrian)) @@ -172,6 +165,10 @@ function CellData.change_domain(a::CellField,strian::AdaptedTriangulation,::Refe end function CellData.change_domain(a::CellField,strian::AdaptedTriangulation,::ReferenceDomain,ttrian::AdaptedTriangulation,::ReferenceDomain) + if strian === ttrian + return a + end + if is_child(strian,ttrian) # fine to coarse b = change_domain(a,strian,ReferenceDomain(),ttrian.trian,ReferenceDomain()) return CellData.similar_cell_field(b,get_data(b),ttrian,ReferenceDomain()) @@ -228,6 +225,7 @@ function change_domain_o2n(f_old,old_trian::Triangulation,new_trian::AdaptedTria end function change_domain_o2n(f_coarse,ctrian::Triangulation{Dc},ftrian::AdaptedTriangulation,glue::AdaptivityGlue{<:RefinementGlue}) where Dc + @notimplementedif num_point_dims(ctrian) != Dc cglue = get_glue(ctrian,Val(Dc)) fglue = get_glue(ftrian,Val(Dc)) @@ -267,6 +265,7 @@ function change_domain_n2o(f_new,new_trian::AdaptedTriangulation,old_trian::Tria end function change_domain_n2o(f_fine,ftrian::AdaptedTriangulation{Dc},ctrian::Triangulation,glue::AdaptivityGlue{<:RefinementGlue}) where Dc + @notimplementedif num_point_dims(ftrian) != Dc cglue = get_glue(ctrian,Val(Dc)) fglue = get_glue(ftrian,Val(Dc)) From 3d16c89781a0ea9dcd42e32e9bcd3c08142dcba9 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 14 Feb 2023 10:12:38 +1100 Subject: [PATCH 6/9] Added tests --- ...inTests.jl => ComplexChangeDomainTests.jl} | 2 +- test/AdaptivityTests/mytests.jl | 89 +++++++++++++++++++ test/AdaptivityTests/runtests.jl | 1 + 3 files changed, 91 insertions(+), 1 deletion(-) rename test/AdaptivityTests/{ExtendingChangeDomainTests.jl => ComplexChangeDomainTests.jl} (98%) create mode 100644 test/AdaptivityTests/mytests.jl diff --git a/test/AdaptivityTests/ExtendingChangeDomainTests.jl b/test/AdaptivityTests/ComplexChangeDomainTests.jl similarity index 98% rename from test/AdaptivityTests/ExtendingChangeDomainTests.jl rename to test/AdaptivityTests/ComplexChangeDomainTests.jl index 242050a26..2e2ff1e6a 100644 --- a/test/AdaptivityTests/ExtendingChangeDomainTests.jl +++ b/test/AdaptivityTests/ComplexChangeDomainTests.jl @@ -1,4 +1,4 @@ -module ExtendingChangeDomainTests +module ComplexChangeDomainTests using Test using Gridap diff --git a/test/AdaptivityTests/mytests.jl b/test/AdaptivityTests/mytests.jl new file mode 100644 index 000000000..d22b82f41 --- /dev/null +++ b/test/AdaptivityTests/mytests.jl @@ -0,0 +1,89 @@ +module MyTests + +using Test +using Gridap +using Gridap.Geometry +using Gridap.CellData +using Gridap.Adaptivity +using Gridap.ReferenceFEs +using FillArrays + +""" +Keep in mind: There is a function _compose_glues() that +we could extend! +""" +sol(x) = sum(x) + +D = 2 +domain = Tuple(repeat([0,1],D)) +cmodel = CartesianDiscreteModel(domain,Tuple(fill(2,D))) +fmodel = refine(cmodel,2) + +Ωc = Triangulation(cmodel) +Ωf = Triangulation(fmodel) + +Λc = Skeleton(cmodel) +Λf = Skeleton(fmodel) + +reffe = ReferenceFE(lagrangian,Float64,1) +Vc = TestFESpace(cmodel,reffe) +Uc = TrialFESpace(Vc) +Vf = TestFESpace(fmodel,reffe) +Uf = TrialFESpace(Vf) + +""" + So this works out of the box for any change between a Dc==Dp triangulation and + any other triangulation for which the change is possible. +""" +v_Ωc = interpolate(sol,Vc) +aux = Triangulation(ReferenceFE{num_cell_dims(Ωc)},get_adapted_model(Λf)) +v_Ωf = change_domain(v_Ωc.plus,aux,ReferenceDomain()) +v_Λf = change_domain(v_Ωf,Λf,ReferenceDomain()) + +v_Ωf = interpolate(sol,Vf) +aux = Triangulation(ReferenceFE{num_cell_dims(Ωf)},get_background_model(Λc)) +v_Ωc = change_domain(v_Ωf.plus,aux,ReferenceDomain()) +v_Λc = change_domain(v_Ωc,Λc,ReferenceDomain()) + + +""" +Attempt at an edge adaptivity glue +""" + +glue = get_adaptivity_glue(fmodel) +n2o_cells = glue.n2o_faces_map[3] +o2n_cells = glue.o2n_faces_map + +ctopo = get_grid_topology(cmodel) +ftopo = get_grid_topology(fmodel) + +c2e_new = get_faces(ftopo,2,1) +c2e_old = get_faces(ctopo,2,1) + +rr = glue.refinement_rules[1] +rr_grid = Adaptivity.get_ref_grid(rr) + +# [n/e][child_id][flid]->clid +const rrule_f_to_c_lid_2D=Vector{Vector{Vector{UInt8}}}( + [ + [[1,1,3,1], [1,2,1,4], [3,1,3,2], [1,4,2,4]], # nodes + [[1,1,3,1], [1,1,1,4], [1,2,3,1], [1,2,1,4]] # edges + ]) + +# [n/e][child_id][flid]->clid_dim +const rrule_f_to_c_dim_2D=Vector{Vector{Vector{UInt8}}}( + [ + [[0,1,1,2], [1,0,2,1], [1,2,0,1], [2,1,1,0]], # nodes + [[1,2,1,2], [1,2,2,1], [2,1,1,2], [2,1,2,1]] # edges + ]) + +nE = num_edges(fmodel) +n2o_edges = -ones(nE) + +iO = 1 +new_cells = o2n_cells[iO] +new_edges = c2e_new[new_cells] +old_edges = c2e_old[iO] + + +n2o_edges[] \ No newline at end of file diff --git a/test/AdaptivityTests/runtests.jl b/test/AdaptivityTests/runtests.jl index 0e3d2495a..588e556e5 100644 --- a/test/AdaptivityTests/runtests.jl +++ b/test/AdaptivityTests/runtests.jl @@ -9,6 +9,7 @@ end @testset "Refinement" begin include("CartesianRefinementTests.jl") + include("ComplexChangeDomainTests.jl") include("EdgeBasedRefinementTests.jl") include("FineToCoarseFieldsTests.jl") end From fffd62fc597c1db501cf30df729e1e0fd7a2cedb Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 14 Feb 2023 10:21:23 +1100 Subject: [PATCH 7/9] Updated NEWS.md --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 15ee7ee6a..92ba35ac2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Implemented `CompositeQuadrature`, a quadrature for a cell that has been refined using a `RefinementRule`. Since PR [#838](https://github.com/gridap/Gridap.jl/pull/838). - Implemented simple refinement strategies for Cartesian discrete models in 2&3D as well as Unstructured discrete models in 2D. The latter is implemented by red-green refinement. Since PR [#838](https://github.com/gridap/Gridap.jl/pull/838). - Added optimization when calling `unique` for a `CompressedArray`. Since PR [#838](https://github.com/gridap/Gridap.jl/pull/838). +- Added support for changing domain between adapted triangulations in cases where the target triangulation is a view, a `BoundaryTriangulation` or a `SkeletonTriangulation`. Since PR [#868](https://github.com/gridap/Gridap.jl/pull/868). ### Fixed - Using broadcasting through in `ODESolver` vector operations. Since PR [#860](https://github.com/gridap/Gridap.jl/pull/860) From d6a53e3173e9ce4981ad78c67c1654975c5d9370 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 21 Feb 2023 14:57:45 +1100 Subject: [PATCH 8/9] adapted move_contributions now supports all triangulation types --- src/Adaptivity/AdaptedTriangulations.jl | 40 ++++++++++++++++--- .../ComplexChangeDomainTests.jl | 10 +++++ 2 files changed, 44 insertions(+), 6 deletions(-) diff --git a/src/Adaptivity/AdaptedTriangulations.jl b/src/Adaptivity/AdaptedTriangulations.jl index f25e16510..b6f234f39 100644 --- a/src/Adaptivity/AdaptedTriangulations.jl +++ b/src/Adaptivity/AdaptedTriangulations.jl @@ -189,12 +189,36 @@ for sdomain in [:ReferenceDomain,:PhysicalDomain] end end -function Geometry.move_contributions(scell_to_val::AbstractArray, strian::AdaptedTriangulation, ttrian::Triangulation) +function Geometry.move_contributions(scell_to_val::AbstractArray, strian::AdaptedTriangulation, ttrian::Triangulation{Dc}) where Dc + # Default Gridap + if get_background_model(strian) === get_background_model(ttrian) + return move_contributions(scell_to_val,strian.trian,ttrian) + end + + # Else smodel = get_adapted_model(strian) @check get_parent(smodel) === get_background_model(ttrian) - tcell_to_val = move_contributions(scell_to_val,get_adaptivity_glue(smodel)) - return tcell_to_val + return move_contributions_f2c(scell_to_val,strian,ttrian) +end + +function move_contributions_f2c(fine_tface_to_val::AbstractArray, ftrian::AdaptedTriangulation, ctrian::Triangulation{Dc}) where Dc + glue = get_adaptivity_glue(ftrian) + fglue = get_glue(ftrian,Val(Dc)) + cglue = get_glue(ctrian,Val(Dc)) + + # Fine Triangulation -> Fine Model + fmodel = get_background_model(ftrian) + cell_ftrian = Triangulation(ReferenceFE{Dc},fmodel) + fine_mface_to_val = move_contributions(fine_tface_to_val,ftrian.trian,cell_ftrian) + + # Fine model -> Coarse Model + coarse_mface_to_val = move_contributions(fine_mface_to_val,glue) + + # Coarse Model -> Coarse Triangulation + coarse_tface_to_val = lazy_map(Reindex(coarse_mface_to_val),cglue.tface_to_mface) + + return coarse_tface_to_val end function Geometry.move_contributions(scell_to_val::AbstractArray, glue::AdaptivityGlue) @@ -225,10 +249,12 @@ function change_domain_o2n(f_old,old_trian::Triangulation,new_trian::AdaptedTria end function change_domain_o2n(f_coarse,ctrian::Triangulation{Dc},ftrian::AdaptedTriangulation,glue::AdaptivityGlue{<:RefinementGlue}) where Dc - @notimplementedif num_point_dims(ctrian) != Dc cglue = get_glue(ctrian,Val(Dc)) fglue = get_glue(ftrian,Val(Dc)) + @notimplementedif num_point_dims(ctrian) != Dc + @notimplementedif isa(fglue,Nothing) + if (num_cells(ctrian) != 0) ### Old Triangulation -> Old Model coarse_tface_to_field = CellData.get_data(f_coarse) @@ -265,9 +291,11 @@ function change_domain_n2o(f_new,new_trian::AdaptedTriangulation,old_trian::Tria end function change_domain_n2o(f_fine,ftrian::AdaptedTriangulation{Dc},ctrian::Triangulation,glue::AdaptivityGlue{<:RefinementGlue}) where Dc - @notimplementedif num_point_dims(ftrian) != Dc - cglue = get_glue(ctrian,Val(Dc)) fglue = get_glue(ftrian,Val(Dc)) + cglue = get_glue(ctrian,Val(Dc)) + + @notimplementedif num_point_dims(ftrian) != Dc + @notimplementedif isa(cglue,Nothing) if (num_cells(ctrian) != 0) ### New Triangulation -> New Model diff --git a/test/AdaptivityTests/ComplexChangeDomainTests.jl b/test/AdaptivityTests/ComplexChangeDomainTests.jl index 2e2ff1e6a..10efe88f6 100644 --- a/test/AdaptivityTests/ComplexChangeDomainTests.jl +++ b/test/AdaptivityTests/ComplexChangeDomainTests.jl @@ -8,6 +8,8 @@ using Gridap.Adaptivity using Gridap.ReferenceFEs using FillArrays +dot_prod(u,v,Ω::Triangulation) = sum(∫(v⋅u)*Measure(Ω,qorder)) +dot_prod(u,v,dΩ::Measure) = sum(∫(v⋅u)*dΩ) norm2(u,Ω::Triangulation) = sum(∫(u⋅u)*Measure(Ω,qorder)) norm2(u,dΩ::Measure) = sum(∫(u⋅u)*dΩ) @@ -53,6 +55,10 @@ v_Ωf_view_view = change_domain(v_Ωc_view,Ωf_view,ReferenceDomain()) @test norm2(v_Ωc_view,Ωf_view) ≈ norm2(v_Ωf_view,Ωf_view) @test norm2(v_Ωf_view,Ωc_view) ≈ norm2(v_Ωc_view,Ωc_view) +# Integrate over Ωc_view by integrating first in Ωf and then moving contributions +dΩ_fc = Measure(Ωc_view,Ωf,qorder) +@test dot_prod(v_Ωf,v_Ωc,dΩ_fc) ≈ dot_prod(v_Ωf,v_Ωc,Ωc_view) + """ BodyFitted -> Boundary """ @@ -68,6 +74,10 @@ v_Γc = change_domain(v_Ωc,Γc,ReferenceDomain()) @test norm2(v_Ωf,Γc) ≈ norm2(v_Γc,Γc) @test norm2(v_Ωc,Γf) ≈ norm2(v_Γf,Γf) +# Integrate in Γf then add contributions for each coarse cell +dΩ_fc = Measure(Ωc,Γf,qorder) +@test dot_prod(v_Ωf,v_Ωc,dΩ_fc) ≈ dot_prod(v_Ωf,v_Ωc,Γf) + """ BodyFitted -> Skeleton From f6928a20934804e1e786bd8b064916b532c02313 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 23 Feb 2023 10:10:45 +1100 Subject: [PATCH 9/9] Cleanup --- test/AdaptivityTests/mytests.jl | 89 --------------------------------- 1 file changed, 89 deletions(-) delete mode 100644 test/AdaptivityTests/mytests.jl diff --git a/test/AdaptivityTests/mytests.jl b/test/AdaptivityTests/mytests.jl deleted file mode 100644 index d22b82f41..000000000 --- a/test/AdaptivityTests/mytests.jl +++ /dev/null @@ -1,89 +0,0 @@ -module MyTests - -using Test -using Gridap -using Gridap.Geometry -using Gridap.CellData -using Gridap.Adaptivity -using Gridap.ReferenceFEs -using FillArrays - -""" -Keep in mind: There is a function _compose_glues() that -we could extend! -""" -sol(x) = sum(x) - -D = 2 -domain = Tuple(repeat([0,1],D)) -cmodel = CartesianDiscreteModel(domain,Tuple(fill(2,D))) -fmodel = refine(cmodel,2) - -Ωc = Triangulation(cmodel) -Ωf = Triangulation(fmodel) - -Λc = Skeleton(cmodel) -Λf = Skeleton(fmodel) - -reffe = ReferenceFE(lagrangian,Float64,1) -Vc = TestFESpace(cmodel,reffe) -Uc = TrialFESpace(Vc) -Vf = TestFESpace(fmodel,reffe) -Uf = TrialFESpace(Vf) - -""" - So this works out of the box for any change between a Dc==Dp triangulation and - any other triangulation for which the change is possible. -""" -v_Ωc = interpolate(sol,Vc) -aux = Triangulation(ReferenceFE{num_cell_dims(Ωc)},get_adapted_model(Λf)) -v_Ωf = change_domain(v_Ωc.plus,aux,ReferenceDomain()) -v_Λf = change_domain(v_Ωf,Λf,ReferenceDomain()) - -v_Ωf = interpolate(sol,Vf) -aux = Triangulation(ReferenceFE{num_cell_dims(Ωf)},get_background_model(Λc)) -v_Ωc = change_domain(v_Ωf.plus,aux,ReferenceDomain()) -v_Λc = change_domain(v_Ωc,Λc,ReferenceDomain()) - - -""" -Attempt at an edge adaptivity glue -""" - -glue = get_adaptivity_glue(fmodel) -n2o_cells = glue.n2o_faces_map[3] -o2n_cells = glue.o2n_faces_map - -ctopo = get_grid_topology(cmodel) -ftopo = get_grid_topology(fmodel) - -c2e_new = get_faces(ftopo,2,1) -c2e_old = get_faces(ctopo,2,1) - -rr = glue.refinement_rules[1] -rr_grid = Adaptivity.get_ref_grid(rr) - -# [n/e][child_id][flid]->clid -const rrule_f_to_c_lid_2D=Vector{Vector{Vector{UInt8}}}( - [ - [[1,1,3,1], [1,2,1,4], [3,1,3,2], [1,4,2,4]], # nodes - [[1,1,3,1], [1,1,1,4], [1,2,3,1], [1,2,1,4]] # edges - ]) - -# [n/e][child_id][flid]->clid_dim -const rrule_f_to_c_dim_2D=Vector{Vector{Vector{UInt8}}}( - [ - [[0,1,1,2], [1,0,2,1], [1,2,0,1], [2,1,1,0]], # nodes - [[1,2,1,2], [1,2,2,1], [2,1,1,2], [2,1,2,1]] # edges - ]) - -nE = num_edges(fmodel) -n2o_edges = -ones(nE) - -iO = 1 -new_cells = o2n_cells[iO] -new_edges = c2e_new[new_cells] -old_edges = c2e_old[iO] - - -n2o_edges[] \ No newline at end of file