diff --git a/NEWS.md b/NEWS.md index 8eae5a9..e78f544 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added support for multiple ghost layers on cartesian models. Since PR[#182](https://github.com/gridap/GridapDistributed.jl/pull/182). +### Fixed + +- Fixed issue [#177](https://github.com/gridap/GridapDistributed.jl/issues/177) and [#170](https://github.com/gridap/GridapDistributed.jl/issues/170). Since PR[#180](https://github.com/gridap/GridapDistributed.jl/pull/180). +- Fixed issue where calling `Boundary(with_ghost, dmodel)` would return the local processor boundaries (which include the faces at the interface between processors) instead of returning the local part of the global boundary. Since PR[#180](https://github.com/gridap/GridapDistributed.jl/pull/180). + ## [0.4.9] - 2025-08-08 ### Added diff --git a/src/Geometry.jl b/src/Geometry.jl index bbda42b..e72e26a 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -38,14 +38,27 @@ Geometry.num_point_dims(::Type{<:DistributedGrid{Dc,Dp}}) where {Dc,Dp} = Dp # This object cannot implement the GridTopology interface in a strict sense """ """ -struct DistributedGridTopology{Dc,Dp,A} <: GridapType +struct DistributedGridTopology{Dc,Dp,A,B} <: GridapType topos::A - function DistributedGridTopology(topos::AbstractArray{<:GridTopology{Dc,Dp}}) where {Dc,Dp} + face_gids::B + function DistributedGridTopology( + topos::AbstractArray{<:GridTopology{Dc,Dp}}, + face_gids::AbstractArray{<:PRange} + ) where {Dc,Dp} A = typeof(topos) - new{Dc,Dp,A}(topos) + B = typeof(face_gids) + new{Dc,Dp,A,B}(topos,face_gids) end end +function DistributedGridTopology( + topos::AbstractArray{<:GridTopology{Dc,Dp}}, cell_gids::PRange +) where {Dc,Dp} + face_gids = Vector{PRange}(undef,Dc+1) + face_gids[Dc+1] = cell_gids + return DistributedGridTopology(topos,face_gids) +end + local_views(a::DistributedGridTopology) = a.topos function Geometry.OrientationStyle( @@ -63,6 +76,39 @@ Geometry.num_cell_dims(::Type{<:DistributedGridTopology{Dc,Dp}}) where {Dc,Dp} = Geometry.num_point_dims(::DistributedGridTopology{Dc,Dp}) where {Dc,Dp} = Dp Geometry.num_point_dims(::Type{<:DistributedGridTopology{Dc,Dp}}) where {Dc,Dp} = Dp +function get_cell_gids(topo::DistributedGridTopology{Dc}) where Dc + topo.face_gids[Dc+1] +end + +function get_face_gids(topo::DistributedGridTopology,dim::Integer) + _setup_face_gids!(topo,dim) + return topo.face_gids[dim+1] +end + +function _setup_face_gids!(topo::DistributedGridTopology{Dc},dim) where {Dc} + Gridap.Helpers.@check 0 <= dim <= Dc + if !isassigned(topo.face_gids,dim+1) + cell_gids = topo.face_gids[Dc+1] + nlfaces = map(local_views(topo)) do topo + num_faces(topo,dim) + end + cell_lfaces = map(local_views(topo)) do topo + get_faces(topo, Dc, dim) + end + topo.face_gids[dim+1] = generate_gids(cell_gids,cell_lfaces,nlfaces) + end +end + +function Geometry.get_isboundary_face(topo::DistributedGridTopology, d::Integer) + face_gids = get_face_gids(topo, d) + is_local_boundary = map(local_views(topo)) do topo + get_isboundary_face(topo,d) + end + t = assemble!(&,PVector(is_local_boundary, partition(face_gids))) + is_global_boundary = partition(fetch(consistent!(fetch(t)))) + return is_global_boundary +end + """ """ struct DistributedFaceLabeling{A<:AbstractArray{<:FaceLabeling}} @@ -72,12 +118,40 @@ end local_views(a::DistributedFaceLabeling) = a.labels function Geometry.add_tag_from_tags!(labels::DistributedFaceLabeling, name, tags) - map(labels.labels) do labels + map(local_views(labels)) do labels add_tag_from_tags!(labels, name, tags) end end -# Dsitributed Discrete models +function Geometry.get_face_mask(labels::DistributedFaceLabeling, tags, d::Integer) + map(local_views(labels)) do labels + get_face_mask(labels, tags, d) + end +end + +function Geometry.FaceLabeling(topo::DistributedGridTopology) + D = num_cell_dims(topo) + labels = map(local_views(topo)) do topo + d_to_ndfaces = [ num_faces(topo,d) for d in 0:D ] + labels = FaceLabeling(d_to_ndfaces) + for d in 0:D + get_face_entity(labels,d) .= 1 # Interior as default + end + add_tag!(labels,"interior",[1]) + add_tag!(labels,"boundary",[2]) + return labels + end + for d in 0:D-1 + dface_to_is_boundary = get_isboundary_face(topo,d) # Global boundary + map(labels,dface_to_is_boundary) do labels, dface_to_is_boundary + dface_to_entity = get_face_entity(labels,d) + dface_to_entity .+= dface_to_is_boundary + end + end + return labels +end + +# Distributed Discrete models # We do not inherit from DiscreteModel on purpose. # This object cannot implement the DiscreteModel interface in a strict sense @@ -131,7 +205,7 @@ function Geometry.get_grid(model::DistributedDiscreteModel) end function Geometry.get_grid_topology(model::DistributedDiscreteModel) - DistributedGridTopology(map(get_grid_topology,local_views(model))) + DistributedGridTopology(map(get_grid_topology,local_views(model)),model.face_gids) end function Geometry.get_face_labeling(model::DistributedDiscreteModel) @@ -458,15 +532,19 @@ end # This object cannot implement the Triangulation interface in a strict sense """ """ -struct DistributedTriangulation{Dc,Dp,A,B} <: GridapType - trians::A - model::B +struct DistributedTriangulation{Dc,Dp,A,B,C} <: GridapType + trians ::A + model ::B + metadata::C function DistributedTriangulation( trians::AbstractArray{<:Triangulation{Dc,Dp}}, - model::DistributedDiscreteModel) where {Dc,Dp} + model::DistributedDiscreteModel; + metadata = nothing + ) where {Dc,Dp} A = typeof(trians) B = typeof(model) - new{Dc,Dp,A,B}(trians,model) + C = typeof(metadata) + new{Dc,Dp,A,B,C}(trians,model,metadata) end end @@ -501,58 +579,84 @@ end # Triangulation constructors -function Geometry.Triangulation( - model::DistributedDiscreteModel;kwargs...) - D=num_cell_dims(model) +function Geometry.Triangulation(model::DistributedDiscreteModel;kwargs...) + D = num_cell_dims(model) Triangulation(no_ghost,ReferenceFE{D},model;kwargs...) end -function Geometry.BoundaryTriangulation( - model::DistributedDiscreteModel;kwargs...) - BoundaryTriangulation(no_ghost,model;kwargs...) +function Geometry.Triangulation(::Type{ReferenceFE{D}},model::DistributedDiscreteModel;kwargs...) where D + Triangulation(no_ghost, ReferenceFE{D}, model; kwargs...) end -function Geometry.BoundaryTriangulation( - trian::DistributedTriangulation;kwargs...) - BoundaryTriangulation(no_ghost,trian;kwargs...) +function Geometry.Triangulation(portion, model::DistributedDiscreteModel;kwargs...) + D = num_cell_dims(model) + Triangulation(portion,ReferenceFE{D},model;kwargs...) end -function Geometry.SkeletonTriangulation( - model::DistributedDiscreteModel;kwargs...) - SkeletonTriangulation(no_ghost,model;kwargs...) +function Geometry.Triangulation( + portion,::Type{ReferenceFE{D}},model::DistributedDiscreteModel;kwargs...) where D + gids = get_face_gids(model,D) + trians = map(local_views(model)) do model + Triangulation(ReferenceFE{D},model;kwargs...) + end + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end -function Geometry.SkeletonTriangulation( - trian::DistributedTriangulation;kwargs...) - SkeletonTriangulation(no_ghost,trian;kwargs...) +function Geometry.BoundaryTriangulation(model::DistributedDiscreteModel,args...;kwargs...) + BoundaryTriangulation(no_ghost,model,args...;kwargs...) end -function Geometry.Triangulation( - portion,::Type{ReferenceFE{Dt}},model::DistributedDiscreteModel{Dm};kwargs...) where {Dt,Dm} - # Generate global ordering for the faces of dimension Dt (if needed) - gids = get_face_gids(model,Dt) - trians = map(local_views(model),partition(gids)) do model, gids - Triangulation(portion,gids,ReferenceFE{Dt},model;kwargs...) +function Geometry.BoundaryTriangulation(trian::DistributedTriangulation;kwargs...) + BoundaryTriangulation(no_ghost,trian;kwargs...) +end + +function Geometry.BoundaryTriangulation( + portion,model::DistributedDiscreteModel;kwargs...) + labels = get_face_labeling(model) + Geometry.BoundaryTriangulation(portion,model,labels;kwargs...) +end + +function Geometry.BoundaryTriangulation( + portion,model::DistributedDiscreteModel,labels::DistributedFaceLabeling;tags=nothing) + Dc = num_cell_dims(model) + if isnothing(tags) + topo = get_grid_topology(model) + face_to_mask = get_isboundary_face(topo,Dc-1) # This is globally consistent + else + face_to_mask = get_face_mask(labels,tags,Dc-1) end - DistributedTriangulation(trians,model) + Geometry.BoundaryTriangulation(portion,model,face_to_mask) end function Geometry.BoundaryTriangulation( - portion,model::DistributedDiscreteModel{Dc};kwargs...) where Dc - gids = get_face_gids(model,Dc) - trians = map(local_views(model),partition(gids)) do model, gids - BoundaryTriangulation(portion,gids,model;kwargs...) + portion,model::DistributedDiscreteModel,face_to_mask::AbstractArray) + Dc = num_cell_dims(model) + gids = get_face_gids(model,Dc) + trians = map(local_views(model),face_to_mask) do model, face_to_mask + BoundaryTriangulation(model,face_to_mask) end - DistributedTriangulation(trians,model) + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) +end + +function Geometry.SkeletonTriangulation(model::DistributedDiscreteModel;kwargs...) + SkeletonTriangulation(no_ghost,model;kwargs...) +end + +function Geometry.SkeletonTriangulation(trian::DistributedTriangulation;kwargs...) + SkeletonTriangulation(no_ghost,trian;kwargs...) end function Geometry.SkeletonTriangulation( - portion,model::DistributedDiscreteModel{Dc};kwargs...) where Dc - gids = get_face_gids(model,Dc) - trians = map(local_views(model),partition(gids)) do model, gids - SkeletonTriangulation(portion,gids,model;kwargs...) + portion,model::DistributedDiscreteModel;kwargs...) + Dc = num_cell_dims(model) + gids = get_face_gids(model,Dc) + trians = map(local_views(model)) do model + SkeletonTriangulation(model;kwargs...) end - DistributedTriangulation(trians,model) + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end # NOTE: The following constructors require adding back the ghost cells: @@ -564,10 +668,11 @@ function Geometry.BoundaryTriangulation( model = get_background_model(trian) gids = get_cell_gids(model) ghosted_trian = add_ghost_cells(trian) - trians = map(local_views(ghosted_trian),partition(gids)) do trian, gids - BoundaryTriangulation(portion,gids,trian;kwargs...) + trians = map(local_views(ghosted_trian)) do trian + BoundaryTriangulation(trian;kwargs...) end - DistributedTriangulation(trians,model) + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end function Geometry.SkeletonTriangulation( @@ -576,32 +681,29 @@ function Geometry.SkeletonTriangulation( model = get_background_model(trian) gids = get_cell_gids(model) ghosted_trian = add_ghost_cells(trian) - trians = map(local_views(ghosted_trian),partition(gids)) do trian, gids - SkeletonTriangulation(portion,gids,trian;kwargs...) + trians = map(local_views(ghosted_trian)) do trian + SkeletonTriangulation(trian;kwargs...) end - DistributedTriangulation(trians,model) + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end -function Geometry.Triangulation( - portion,gids::AbstractLocalIndices, args...;kwargs...) +function Geometry.Triangulation(portion,gids::AbstractLocalIndices, args...;kwargs...) trian = Triangulation(args...;kwargs...) filter_cells_when_needed(portion,gids,trian) end -function Geometry.BoundaryTriangulation( - portion,gids::AbstractLocalIndices,args...;kwargs...) +function Geometry.BoundaryTriangulation(portion,gids::AbstractLocalIndices,args...;kwargs...) trian = BoundaryTriangulation(args...;kwargs...) filter_cells_when_needed(portion,gids,trian) end -function Geometry.SkeletonTriangulation( - portion,gids::AbstractLocalIndices,args...;kwargs...) +function Geometry.SkeletonTriangulation(portion,gids::AbstractLocalIndices,args...;kwargs...) trian = SkeletonTriangulation(args...;kwargs...) filter_cells_when_needed(portion,gids,trian) end -function Geometry.InterfaceTriangulation( - portion,gids::AbstractLocalIndices,args...;kwargs...) +function Geometry.InterfaceTriangulation(portion,gids::AbstractLocalIndices,args...;kwargs...) trian = InterfaceTriangulation(args...;kwargs...) filter_cells_when_needed(portion,gids,trian) end @@ -612,53 +714,29 @@ function Geometry.InterfaceTriangulation(a::DistributedTriangulation,b::Distribu DistributedTriangulation(trians,a.model) end -function Geometry.Triangulation( - portion, model::DistributedDiscreteModel;kwargs...) - D = num_cell_dims(model) - Triangulation(portion,ReferenceFE{D},model;kwargs...) -end +# Filtering cells -function Geometry.Triangulation( - ::Type{ReferenceFE{D}}, model::DistributedDiscreteModel;kwargs...) where D - Triangulation(no_ghost, ReferenceFE{D}, model; kwargs...) +@inline function filter_cells_when_needed( + portion::Union{WithGhost,FullyAssembledRows},cell_gids,trian) + return trian end -function filter_cells_when_needed( - portion::WithGhost, - cell_gids::AbstractLocalIndices, - trian::Triangulation) - - trian +@inline function filter_cells_when_needed( + portion::Union{NoGhost,SubAssembledRows},cell_gids,trian) + return remove_ghost_cells(trian,cell_gids) end -function filter_cells_when_needed( - portion::NoGhost, - cell_gids::AbstractLocalIndices, - trian::Triangulation) +# Removing ghost cells - remove_ghost_cells(trian,cell_gids) -end - -function filter_cells_when_needed( - portion::FullyAssembledRows, - cell_gids::AbstractLocalIndices, - trian::Triangulation) - - trian -end - -function filter_cells_when_needed( - portion::SubAssembledRows, - cell_gids::AbstractLocalIndices, - trian::Triangulation) - - remove_ghost_cells(trian,cell_gids) +struct RemoveGhostsMetadata{A} + parents::A end function remove_ghost_cells(trian::DistributedTriangulation,gids) trians = map(remove_ghost_cells,local_views(trian),partition(gids)) model = get_background_model(trian) - return DistributedTriangulation(trians,model) + metadata = RemoveGhostsMetadata(local_views(trian)) + return DistributedTriangulation(trians,model;metadata) end function remove_ghost_cells(trian::Triangulation,gids) @@ -678,10 +756,8 @@ function remove_ghost_cells( remove_ghost_cells(glue,trian,gids) end -function remove_ghost_cells( - trian::AdaptedTriangulation{Dc,Dp,<:Union{SkeletonTriangulation,BoundaryTriangulation}},gids -) where {Dc,Dp} - remove_ghost_cells(trian.trian,gids) +function remove_ghost_cells(trian::AdaptedTriangulation,gids) + AdaptedTriangulation(remove_ghost_cells(trian.trian,gids),trian.adapted_model) end function remove_ghost_cells(glue::FaceToFaceGlue,trian,gids) @@ -724,29 +800,30 @@ function _find_owned_skeleton_facets(glue,gids) end tface_to_part[tface] = part end - findall(part->part==part_id(gids),tface_to_part) + return findall(isequal(part_id(gids)),tface_to_part) end +# Adding ghost cells + function add_ghost_cells(dtrian::DistributedTriangulation) dmodel = get_background_model(dtrian) add_ghost_cells(dmodel,dtrian) end -function _covers_all_faces( - dmodel::DistributedDiscreteModel{Dm}, - dtrian::DistributedTriangulation{Dt} +function add_ghost_cells(dmodel::DistributedDiscreteModel,dtrian::DistributedTriangulation) + add_ghost_cells(dtrian.metadata,dmodel,dtrian) +end + +# We already have the parents saved up +function add_ghost_cells( + metadata::RemoveGhostsMetadata, dmodel::DistributedDiscreteModel{Dm}, dtrian::DistributedTriangulation{Dt} ) where {Dm,Dt} - covers_all_faces = map(local_views(dmodel),local_views(dtrian)) do model, trian - glue = get_glue(trian,Val(Dt)) - @assert isa(glue,FaceToFaceGlue) - isa(glue.tface_to_mface,IdentityVector) - end - reduce(&,covers_all_faces,init=true) + DistributedTriangulation(metadata.parents,dmodel) end +# We have to reconstruct the ghosted triangulation function add_ghost_cells( - dmodel::DistributedDiscreteModel{Dm}, - dtrian::DistributedTriangulation{Dt} + metadata, dmodel::DistributedDiscreteModel{Dm}, dtrian::DistributedTriangulation{Dt} ) where {Dm,Dt} tface_to_mface = map(local_views(dtrian)) do trian @@ -783,6 +860,20 @@ function add_ghost_cells( return DistributedTriangulation(trians,dmodel) end +function _covers_all_faces( + dmodel::DistributedDiscreteModel{Dm}, + dtrian::DistributedTriangulation{Dt} +) where {Dm,Dt} + covers_all_faces = map(local_views(dmodel),local_views(dtrian)) do model, trian + glue = get_glue(trian,Val(Dt)) + @assert isa(glue,FaceToFaceGlue) + isa(glue.tface_to_mface,IdentityVector) + end + reduce(&,covers_all_faces,init=true) +end + +# Triangulation gids + function generate_cell_gids(dtrian::DistributedTriangulation) dmodel = get_background_model(dtrian) generate_cell_gids(dmodel,dtrian)