From d510d2b74da7f9fb628e1821ab9547b0a7e8c62d Mon Sep 17 00:00:00 2001 From: Pere Antoni Martorell Date: Tue, 16 Apr 2024 11:08:37 +0200 Subject: [PATCH] [distributed] missing cut_facets interfaces --- src/Distributed/Distributed.jl | 10 + src/Distributed/DistributedDiscretizations.jl | 106 ++++++++- src/Distributed/DistributedQuadratures.jl | 18 ++ .../MomentFittedQuadratures.jl | 32 ++- test/DistributedTests/DistributedDevs.jl | 214 ------------------ test/DistributedTests/PoissonTests.jl | 4 +- 6 files changed, 156 insertions(+), 228 deletions(-) create mode 100644 src/Distributed/DistributedQuadratures.jl delete mode 100644 test/DistributedTests/DistributedDevs.jl diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl index 54d7669..aadb4d5 100644 --- a/src/Distributed/Distributed.jl +++ b/src/Distributed/Distributed.jl @@ -6,6 +6,7 @@ using PartitionedArrays using FillArrays using Gridap.Arrays +using Gridap.CellData using Gridap.Geometry using Gridap.Helpers using Gridap.ReferenceFEs @@ -18,13 +19,17 @@ using GridapEmbedded.Interfaces: ActiveInOrOut using GridapEmbedded.Interfaces: SubFacetTriangulation using GridapEmbedded.Interfaces: SubCellData using GridapEmbedded.Interfaces: SubFacetData +using GridapEmbedded.Interfaces: AbstractEmbeddedDiscretization using GridapEmbedded.AgFEM: _touch_aggregated_cells! using GridapEmbedded.AgFEM: AggregateCutCellsByThreshold +using GridapEmbedded.MomentFittedQuadratures: MomentFitted using Gridap.Geometry: AppendedTriangulation +using Gridap.Geometry: get_face_to_parent_face using GridapDistributed: DistributedDiscreteModel using GridapDistributed: DistributedTriangulation using GridapDistributed: DistributedFESpace using GridapDistributed: DistributedSingleFieldFESpace +using GridapDistributed: DistributedMeasure using GridapDistributed: add_ghost_cells using GridapDistributed: generate_gids using GridapDistributed: generate_cell_gids @@ -33,11 +38,14 @@ using GridapDistributed: _find_vector_type import GridapEmbedded.AgFEM: aggregate import GridapEmbedded.AgFEM: AgFEMSpace import GridapEmbedded.Interfaces: cut +import GridapEmbedded.Interfaces: cut_facets import GridapEmbedded.Interfaces: EmbeddedBoundary import GridapEmbedded.Interfaces: compute_bgfacet_to_inoutcut import GridapEmbedded.Interfaces: compute_bgcell_to_inoutcut import GridapEmbedded.CSG: get_geometry import Gridap.Geometry: Triangulation +import Gridap.Geometry: SkeletonTriangulation +import Gridap.Geometry: BoundaryTriangulation import Gridap.Geometry: get_background_model import GridapDistributed: local_views import GridapDistributed: remove_ghost_cells @@ -46,4 +54,6 @@ include("DistributedDiscretizations.jl") include("DistributedAgFEM.jl") +include("DistributedQuadratures.jl") + end # module diff --git a/src/Distributed/DistributedDiscretizations.jl b/src/Distributed/DistributedDiscretizations.jl index e2a20a5..850eb18 100644 --- a/src/Distributed/DistributedDiscretizations.jl +++ b/src/Distributed/DistributedDiscretizations.jl @@ -1,13 +1,13 @@ -struct DistributedEmbeddedDiscretization{Dp,T,A,B} <: GridapType +struct DistributedEmbeddedDiscretization{A,B} <: GridapType discretizations::A model::B function DistributedEmbeddedDiscretization( - d::AbstractArray{<:EmbeddedDiscretization{Dp,T}}, - model::DistributedDiscreteModel) where {Dp,T} + d::AbstractArray{<:AbstractEmbeddedDiscretization}, + model::DistributedDiscreteModel) A = typeof(d) B = typeof(model) - new{Dp,T,A,B}(d,model) + new{A,B}(d,model) end end @@ -31,8 +31,29 @@ function cut(cutter::Cutter,bgmodel::DistributedDiscreteModel,args...) cutgeo = cut(cutter,ownmodel,args...) change_bgmodel(cutgeo,bgmodel,own_to_local(gids)) end - ls_to_bgcell_to_inoutcut = map(c->c.ls_to_bgcell_to_inoutcut,cuts) - _consistent!(ls_to_bgcell_to_inoutcut,gids) + consistent_bgcell_to_inoutcut!(cuts,gids) + DistributedEmbeddedDiscretization(cuts,bgmodel) +end + +function cut_facets(bgmodel::DistributedDiscreteModel,args...) + cut_facets(LevelSetCutter(),bgmodel,args...) +end + +function cut_facets(cutter::Cutter,bgmodel::DistributedDiscreteModel,args...) + D = map(num_dims,local_views(bgmodel)) |> PartitionedArrays.getany + cell_gids = get_cell_gids(bgmodel) + facet_gids = get_face_gids(bgmodel,D-1) + cuts = map( + local_views(bgmodel), + local_views(cell_gids), + local_views(facet_gids)) do bgmodel,cell_gids,facet_gids + ownmodel = remove_ghost_cells(bgmodel,cell_gids) + facet_to_pfacet = get_face_to_parent_face(ownmodel,D-1) + cutfacets = cut_facets(cutter,ownmodel,args...) + cutfacets = change_bgmodel(cutfacets,bgmodel,facet_to_pfacet) + remove_ghost_subfacets(cutfacets,facet_gids) + end + consistent_bgfacet_to_inoutcut!(cuts,facet_gids) DistributedEmbeddedDiscretization(cuts,bgmodel) end @@ -54,6 +75,16 @@ function EmbeddedBoundary(cutgeo::DistributedEmbeddedDiscretization,args...) remove_ghost_cells(trian) end +function SkeletonTriangulation(cutgeo::DistributedEmbeddedDiscretization,args...) + trian = distributed_embedded_triangulation(SkeletonTriangulation,cutgeo,args...) + remove_ghost_cells(trian) +end + +function BoundaryTriangulation(cutgeo::DistributedEmbeddedDiscretization,args...) + trian = distributed_embedded_triangulation(BoundaryTriangulation,cutgeo,args...) + remove_ghost_cells(trian) +end + function distributed_embedded_triangulation( T, cutgeo::DistributedEmbeddedDiscretization, @@ -139,6 +170,30 @@ function remove_ghost_cells(model::DiscreteModel,gids::AbstractLocalIndices) DiscreteModelPortion(model,own_to_local(gids)) end +function consistent_bgcell_to_inoutcut!( + cuts::AbstractArray{<:AbstractEmbeddedDiscretization}, + gids::PRange) + + ls_to_bgcell_to_inoutcut = map(get_ls_to_bgcell_to_inoutcut,cuts) + _consistent!(ls_to_bgcell_to_inoutcut,gids) +end + +function get_ls_to_bgcell_to_inoutcut(cut::EmbeddedDiscretization) + cut.ls_to_bgcell_to_inoutcut +end + +function consistent_bgfacet_to_inoutcut!( + cuts::AbstractArray{<:AbstractEmbeddedDiscretization}, + gids::PRange) + + ls_to_bgfacet_to_inoutcut = map(get_ls_to_bgfacet_to_inoutcut,cuts) + _consistent!(ls_to_bgfacet_to_inoutcut,gids) +end + +function get_ls_to_bgfacet_to_inoutcut(cut::EmbeddedFacetDiscretization) + cut.ls_to_facet_to_inoutcut +end + function _consistent!( p_to_i_to_a::AbstractArray{<:Vector{<:Vector}}, prange::PRange) @@ -207,6 +262,28 @@ function change_bgmodel( cut.geo) end +function change_bgmodel( + cut::EmbeddedFacetDiscretization, + newmodel::DiscreteModel, + facet_to_newfacet=1:num_facets(get_background_model(cut))) + + nfacets = num_facets(newmodel) + + ls_to_bgf_to_ioc = map(cut.ls_to_facet_to_inoutcut) do bgf_to_ioc + new_bgf_to_ioc = Vector{Int8}(undef,nfacets) + new_bgf_to_ioc[facet_to_newfacet] = bgf_to_ioc + new_bgf_to_ioc + end + subfacets = change_bgmodel(cut.subfacets,facet_to_newfacet) + EmbeddedFacetDiscretization( + newmodel, + ls_to_bgf_to_ioc, + subfacets, + cut.ls_to_subfacet_to_inout, + cut.oid_to_ls, + cut.geo) +end + function change_bgmodel(cells::SubCellData,cell_to_newcell) cell_to_bgcell = lazy_map(Reindex(cell_to_newcell),cells.cell_to_bgcell) SubCellData( @@ -225,3 +302,20 @@ function change_bgmodel(facets::SubFacetData,cell_to_newcell) facets.point_to_coords, facets.point_to_rcoords) end + +function remove_ghost_subfacets(cut::EmbeddedFacetDiscretization,facet_gids) + bgfacet_mask = map(!iszero,local_to_owner(facet_gids)) + subfacet_mask = map(Reindex(bgfacet_mask),cut.subfacets.cell_to_bgcell) + new_subfacets = findall(subfacet_mask) + subfacets = SubCellData(cut.subfacets,new_subfacets) + ls_to_subfacet_to_inout = map(cut.ls_to_subfacet_to_inout) do sf_to_io + map(Reindex(sf_to_io),new_subfacets) + end + EmbeddedFacetDiscretization( + cut.bgmodel, + cut.ls_to_facet_to_inoutcut, + subfacets, + ls_to_subfacet_to_inout, + cut.oid_to_ls, + cut.geo) +end diff --git a/src/Distributed/DistributedQuadratures.jl b/src/Distributed/DistributedQuadratures.jl new file mode 100644 index 0000000..06a40a4 --- /dev/null +++ b/src/Distributed/DistributedQuadratures.jl @@ -0,0 +1,18 @@ + +function CellData.Measure( + t::DistributedTriangulation, + quad::Tuple{MomentFitted,Vararg}; + kwargs...) + @notimplemented + name, _args, _kwargs = quad + cut,cutfacets,_args... = _args + t = remove_ghost_cells(t) + measures = map( + local_views(t), + local_views(cut), + local_views(cutfacets)) do trian,cut,cutfacets + quad = name, (cut,cutfacets,_args...), _kwargs + Measure(trian,quad;kwargs...) + end + DistributedMeasure(measures) +end diff --git a/src/MomentFittedQuadratures/MomentFittedQuadratures.jl b/src/MomentFittedQuadratures/MomentFittedQuadratures.jl index 5c6f114..2102263 100644 --- a/src/MomentFittedQuadratures/MomentFittedQuadratures.jl +++ b/src/MomentFittedQuadratures/MomentFittedQuadratures.jl @@ -43,15 +43,36 @@ function Quadrature(trian::Grid, Quadrature(trian,momentfitted,cut,geo,degree;in_or_out=in_or_out) end +function Quadrature(trian::Grid, + ::MomentFitted, + cut::AbstractEmbeddedDiscretization, + cutf::AbstractEmbeddedDiscretization, + degree::Int; + in_or_out=IN) + geo = get_geometry(cut) + Quadrature(trian,momentfitted,cut,cutf,geo,degree;in_or_out=in_or_out) +end + +function Quadrature(trian::Grid, + ::MomentFitted, + cut::AbstractEmbeddedDiscretization, + geo::CSG.Geometry, + degree::Int; + in_or_out=IN) + cutf = cut_facets(cut,geo) + Quadrature(trian,momentfitted,cut,cutf,geo,degree;in_or_out=in_or_out) +end + function Quadrature(active_mesh::Grid, ::MomentFitted, cut::AbstractEmbeddedDiscretization, + cutf::AbstractEmbeddedDiscretization, geo::CSG.Geometry, degree::Int; in_or_out=IN) acell_to_point_vals, acell_to_weight_vals = # - compute_lag_moments_from_leg(cut,geo,in_or_out,degree) + compute_lag_moments_from_leg(cut,cutf,geo,in_or_out,degree) acell_to_weight_vals = collect(get_array(acell_to_weight_vals)) D = num_dims(active_mesh) @@ -131,6 +152,7 @@ function legendreToMonomial(n::Int,d::Int) end function compute_lag_moments_from_leg(cut::AbstractEmbeddedDiscretization, + cutf::AbstractEmbeddedDiscretization, geo::CSG.Geometry, in_or_out, degree::Int) @@ -138,7 +160,7 @@ function compute_lag_moments_from_leg(cut::AbstractEmbeddedDiscretization, T = eltype(eltype(get_node_coordinates(cut_trian))) D = num_dims(cut_trian) b = MonomialBasis{D}(T,degree) - mon_contribs = compute_monomial_domain_contribution(cut,geo,in_or_out,b,degree) + mon_contribs = compute_monomial_domain_contribution(cut,cutf,geo,in_or_out,b,degree) mon_moments = compute_monomial_cut_cell_moments(cut_trian,mon_contribs,b) mon_to_leg = Fill(legendreToMonomial(degree,D),num_cells(cut_trian)) leg_moments = lazy_map(*,mon_to_leg,mon_moments) @@ -149,14 +171,16 @@ function compute_lag_moments_from_leg(cut::AbstractEmbeddedDiscretization, end function compute_monomial_domain_contribution(cut::AbstractEmbeddedDiscretization, + cutf::AbstractEmbeddedDiscretization, in_or_out, b::MonomialBasis, deg::Int) geo = get_geometry(cut) - compute_monomial_domain_contribution(cut,geo,in_or_out,b,deg) + compute_monomial_domain_contribution(cut,cutf,geo,in_or_out,b,deg) end function compute_monomial_domain_contribution(cut::AbstractEmbeddedDiscretization, + cutf::AbstractEmbeddedDiscretization, geo::CSG.Geometry, in_or_out::Integer, b::MonomialBasis, @@ -164,7 +188,6 @@ function compute_monomial_domain_contribution(cut::AbstractEmbeddedDiscretizatio cut_io = CutInOrOut(in_or_out) dir_Γᵉ = (-1)^(in_or_out==OUT) - cutf = cut_facets(cut,geo) # Embedded facets Γᵉ = EmbeddedBoundary(cut,geo) # Interior fitted cut facets @@ -354,4 +377,3 @@ function _get_terms_degrees(c::CartesianIndex) end end #module - diff --git a/test/DistributedTests/DistributedDevs.jl b/test/DistributedTests/DistributedDevs.jl deleted file mode 100644 index 33f686e..0000000 --- a/test/DistributedTests/DistributedDevs.jl +++ /dev/null @@ -1,214 +0,0 @@ -module DistributedDevs - -using Gridap -using GridapEmbedded -using GridapDistributed -using PartitionedArrays -using Test - -using GridapEmbedded.Distributed: distributed_aggregate -using GridapEmbedded.Distributed: has_remote_aggregation - -distribute = PartitionedArrays.DebugArray - -np = (2,1) - -ranks = distribute(LinearIndices((prod(np),))) - -L = 1 -p0 = Point(0.0,0.0) -pmin = p0-L/2 -pmax = p0+L/2 - -R = 0.3 -x = VectorValue(1.0,0.0) - -d1 = VectorValue(0.4,0.0) -d2 = VectorValue(0.0,0.05) -geo = disk(R,x0=p0) -# geo2 = quadrilateral(;x0=p0-d2/2,d1=d1,d2=d2) -# geo = union(geo1,geo2) - -n = 9 -mesh_partition = (n,n) -bgmodel = CartesianDiscreteModel(ranks,np,pmin,pmax,mesh_partition) - -cutgeo = cut(bgmodel,geo) - -bgf_to_ioc = compute_bgfacet_to_inoutcut(bgmodel,geo) - -Ω = Triangulation(cutgeo) - -writevtk(Ω,"trian") - -strategy = AggregateCutCellsByThreshold(1.0) -aggregates,aggregate_owner,aggregate_neig = distributed_aggregate( - strategy,cutgeo,geo,IN) - - -gids = get_cell_gids(bgmodel) -# TODO: do not use owner for this query (has ghost aggregation?) -local_aggretation = map(aggregate_owner,own_to_local(gids),own_to_owner(gids)) do owns,o_to_l,o - owners = map(Reindex(owns),o_to_l) - all(lazy_map(==,owners,o)) -end -has_local_aggretation = reduction(&,local_aggretation,destination=:all) |> PartitionedArrays.getany - - -@test !has_remote_aggregation(bgmodel,aggregates) - - -@test !has_local_aggretation -# @test !has_remote_aggregation - -cut_in = map(aggregates) do agg - findall(!iszero,agg) -end - -cut_owner = map(aggregate_owner,cut_in) do owners,cut_in - owners[cut_in] -end - -aggregates -gids = get_cell_gids(bgmodel) - -laggregates = map(aggregates,global_to_local(gids)) do agg,g_to_l - map(agg) do a - if iszero(a) - a - else - g_to_l[a] - end - end -end - -oaggregates = map(aggregates,own_to_local(gids)) do agg,o_to_l - map(Reindex(agg),o_to_l) -end - -oaggregate_owner = map(aggregate_owner,own_to_local(gids)) do agg,o_to_l - map(Reindex(agg),o_to_l) -end - - -Ωbg = Triangulation(bgmodel) -Ω = Triangulation(cutgeo) -Ωin = Triangulation(cutgeo,IN) -Γ = EmbeddedBoundary(cutgeo) -Ω_act = Triangulation(cutgeo,ACTIVE) - - - -n_Γ = get_normal_vector(Γ) - -order = 1 -degree = 2*order -dΩ = Measure(Ω,degree) -dΓ = Measure(Γ,degree) - -reffe = ReferenceFE(lagrangian,Float64,order) - -Vstd = FESpace(Ω_act,reffe) - - -map(num_cells,local_views(Ω_act)) - - -map(local_views(Ω_act),ranks) do trian,p - writevtk(trian,"trian_act_$p") -end - -V = AgFEMSpace(bgmodel,Vstd,laggregates) -# V = Vstd -U = TrialFESpace(V) - -h = 1.0/n -γd = 10.0 - -u(x) = x[1] - x[2] -f(x) = -Δ(u)(x) -ud(x) = u(x) - - -a(u,v) = - ∫( ∇(v)⋅∇(u) ) * dΩ + - ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ - -l(v) = - ∫( v*f ) * dΩ + - ∫( (γd/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ - -op = AffineFEOperator(a,l,U,V) -uh = solve(op) - -e = u - uh - -l2(u) = sqrt(sum( ∫( u*u )*dΩ )) -h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ )) - -el2 = l2(e) -eh1 = h1(e) -ul2 = l2(uh) -uh1 = h1(uh) - - -writevtk(Ωin,"trian_in") -writevtk(Γ,"bnd") -writevtk(Ωbg,"bgtrian",celldata= - ["aggregate"=>oaggregates, - "aggregate_owner"=>oaggregate_owner]) - -writevtk(Ω,"trian",cellfields=["uh"=>uh]) - - - - -using GridapEmbedded.Distributed: add_remote_aggregates - -np = (3,) - -ranks = distribute(LinearIndices((prod(np),))) -m = CartesianDiscreteModel(ranks,np,(0, 1),(9,)) - - -aggregates = DebugArray([[1, 9, 8, 5], [3, 4, 9, 6, 7], [6, 7, 3, 9]]) -aggregate_owner = DebugArray([[1, 3, 3, 2], [1, 2, 3, 2, 3], [2, 3, 1, 3]]) - -am = add_remote_aggregates(m,aggregates,aggregate_owner) - -gids = partition(get_cell_gids(am)) -test_gids = DebugArray([[1,2,3,4,5,8,9],[3,4,5,6,7,9],[6,7,8,9,3]]) -map(gids,test_gids) do gids,_gids - @test all(gids .==_gids) -end - -map(local_views(am),ranks) do m,p - writevtk(Triangulation(m),"m$p") -end - - -# import GridapEmbedded.Distributed: change_bgmodel -# using GridapEmbedded.Distributed: DistributedEmbeddedDiscretization -# using GridapDistributed: DistributedDiscreteModel - -# TODO: -# - print aggregates [x] -# - move to src [x] -# - test aggregates with several geometries [x] -# - reconstruct paths [x] -# - add remote roots to model [x] -# - migrate model [x] -# - agfem fe space with remotes - -# TODO: parallel aggregation (parallel aggfem article) -# 0. exchange measures and aggregates through gid [x] -# 1. root in ghost layer [x] -# 2. root in neighbors [x] -# 3. root in neighbors of neighbors [x] - -# TODO: -# - test possion eq with root in ghost layer (simplify this geometry) [x] -# - test graph reconstruction in 1D arrays [x] -# - add remote roots to model [x] - -end diff --git a/test/DistributedTests/PoissonTests.jl b/test/DistributedTests/PoissonTests.jl index 1dd27bf..4c736c3 100644 --- a/test/DistributedTests/PoissonTests.jl +++ b/test/DistributedTests/PoissonTests.jl @@ -20,14 +20,11 @@ function main(distribute,parts; f(x) = -Δ(u)(x) ud(x) = u(x) - - geometries = Dict( :circle => circle_geometry, :remotes => remotes_geometry, ) - bgmodel,geo = geometries[geometry](ranks,parts,cells) D = 2 @@ -36,6 +33,7 @@ function main(distribute,parts; h = meas^(1/D) cutgeo = cut(bgmodel,geo) + cutgeo_facets = cut_facets(bgmodel,geo) strategy = AggregateCutCellsByThreshold(threshold) bgmodel,cutgeo,aggregates = aggregate(strategy,cutgeo)