Constant fespaces

NEWS.md

@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 - Added support for automatic differentiation with ForwardDiff. Since PR[#167](https://github.com/gridap/GridapDistributed.jl/pull/167).
+- Added ConstantFESpaces. Since PR[#166](https://github.com/gridap/GridapDistributed.jl/pull/166).
 
 ## [0.4.5] 2024-10-08
 

src/Adaptivity.jl

@@ -3,6 +3,12 @@
 
 const DistributedAdaptedDiscreteModel{Dc,Dp} = GenericDistributedDiscreteModel{Dc,Dp,<:AbstractArray{<:AdaptedDiscreteModel{Dc,Dp}}}
 
+struct DistributedAdaptedDiscreteModelCache{A,B,C}
+  model_metadata::A
+  parent_metadata::B
+  parent_gids::C
+end
+
 function DistributedAdaptedDiscreteModel(
   model  :: DistributedDiscreteModel,
   parent :: DistributedDiscreteModel,
@@ -12,29 +18,36 @@ function DistributedAdaptedDiscreteModel(
     AdaptedDiscreteModel(model,parent,glue)
   end
   gids = get_cell_gids(model)
-  metadata = hasproperty(model,:metadata) ? model.metadata : nothing
+  metadata = DistributedAdaptedDiscreteModelCache(
+    model.metadata,parent.metadata,get_cell_gids(parent)
+  )
   return GenericDistributedDiscreteModel(models,gids;metadata)
 end
 
 function Adaptivity.get_model(model::DistributedAdaptedDiscreteModel)
   GenericDistributedDiscreteModel(
     map(get_model,local_views(model)),
     get_cell_gids(model);
-    metadata = hasproperty(model,:metadata) ? model.metadata : nothing
+    metadata = model.metadata.model_metadata
   )
 end
 
 function Adaptivity.get_parent(model::DistributedAdaptedDiscreteModel)
-  msg = " Error: Cannot get global parent model. \n 
-          We do not keep the global ids of the parent model within the children.\n
-          You can extract the local parents with map(get_parent,local_views(model))"
-  @notimplemented msg
+  GenericDistributedDiscreteModel(
+    map(get_parent,local_views(model)),
+    model.metadata.parent_gids;
+    metadata = model.metadata.parent_metadata
+  )
 end
 
 function Adaptivity.get_adaptivity_glue(model::DistributedAdaptedDiscreteModel)
   return map(Adaptivity.get_adaptivity_glue,local_views(model))
 end
 
+function Adaptivity.is_child(m1::DistributedDiscreteModel,m2::DistributedDiscreteModel)
+  reduce(&,map(Adaptivity.is_child,local_views(m1),local_views(m2)))
+end
+
 function Adaptivity.refine(
   cmodel::DistributedAdaptedDiscreteModel{Dc},args...;kwargs...
 ) where Dc
@@ -48,7 +61,7 @@ function Adaptivity.refine(
   fmodel = GenericDistributedDiscreteModel(
     map(get_model,local_views(_fmodel)),
     get_cell_gids(_fmodel);
-    metadata=_fmodel.metadata
+    metadata=_fmodel.metadata.model_metadata
   )
   glues = get_adaptivity_glue(_fmodel)
   return DistributedAdaptedDiscreteModel(fmodel,cmodel,glues)
@@ -132,11 +145,7 @@ end
 function redistribute(model::DistributedAdaptedDiscreteModel,args...;kwargs...)
   # Local cmodels are AdaptedDiscreteModels. To correctly dispatch, we need to
   # extract the underlying models, then redistribute.
-  _model = GenericDistributedDiscreteModel(
-    map(get_model,local_views(model)),
-    get_cell_gids(model);
-    metadata=model.metadata
-  )
+  _model = get_model(model)
   return redistribute(_model,args...;kwargs...)
 end
 
@@ -507,7 +516,7 @@ end
 # Cartesian Model uniform refinement
 
 function Adaptivity.refine(
-  cmodel::DistributedDiscreteModel{Dc},
+  cmodel::DistributedCartesianDiscreteModel{Dc},
   refs::Integer = 2
 ) where Dc
   Adaptivity.refine(cmodel,Tuple(fill(refs,Dc)))
@@ -592,7 +601,9 @@ function Adaptivity.refine(
   end
 
   fgids = get_cell_gids(fmodel)
-  metadata = fmodel.metadata
+  metadata = DistributedAdaptedDiscreteModelCache(
+    fmodel.metadata,cmodel.metadata,get_cell_gids(cmodel)
+  )
   return GenericDistributedDiscreteModel(fmodels,fgids;metadata)
 end
 
@@ -909,7 +920,10 @@ function Adaptivity.refine(
 ) where Dc
   fmodels, f_own_to_local = refine_local_models(cmodel,args...;kwargs...)
   fgids = refine_cell_gids(cmodel,fmodels,f_own_to_local)
-  return GenericDistributedDiscreteModel(fmodels,fgids)
+  metadata = DistributedAdaptedDiscreteModelCache(
+    nothing,cmodel.metadata,get_cell_gids(cmodel)
+  )
+  return GenericDistributedDiscreteModel(fmodels,fgids;metadata)
 end
 
 """
@@ -978,7 +992,9 @@ function refine_local_models(
 
   # Filter out local models
   filtered_fmodels = map(fmodels,f_own_or_ghost_ids) do fmodel,f_own_or_ghost_ids
-    model = DiscreteModelPortion(get_model(fmodel),f_own_or_ghost_ids).model
+    model = UnstructuredDiscreteModel( # Necessary to keep the same type
+      DiscreteModelPortion(get_model(fmodel),f_own_or_ghost_ids)
+    )
     parent = get_parent(fmodel)
 
     _glue = get_adaptivity_glue(fmodel)

src/CellData.jl

@@ -362,10 +362,9 @@ struct DistributedInterpolable{Tx,Ty,A} <: Function
   interps::A
   function DistributedInterpolable(interps::AbstractArray{<:Interpolable})
     Tx,Ty = map(interps) do I
-      fi = I.uh
-      trian = get_triangulation(fi)
+      trian = get_triangulation(I.uh)
       x = mean(testitem(get_cell_coordinates(trian)))
-      return typeof(x), return_type(fi,x)
+      return typeof(x), return_type(I,x)
     end |> tuple_of_arrays
     Tx = getany(Tx)
     Ty = getany(Ty)
@@ -418,7 +417,7 @@ function Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVec
       end
       if inside
         ids[k] = i
-        vals[k] = yi
+        vals[k] = copy(yi)
         k += 1
       end
     end

src/FESpaces.jl

@@ -869,3 +869,51 @@ function _compute_new_distributed_fixedval(
   c = reduce(+,c_i,init=zero(eltype(c_i)))
   return c
 end
+
+"""
+    ConstantFESpace(
+      model::DistributedDiscreteModel; 
+      constraint_type=:global, 
+      kwargs...
+    )
+
+Distributed equivalent to `ConstantFESpace(model;kwargs...)`.
+
+With `constraint_type=:global`, a single dof is shared by all processors.
+This creates a global constraint, which is NOT scalable in parallel. Use at your own peril. 
+
+With `constraint_type=:local`, a single dof is owned by each processor and shared with no one else.
+This space is locally-constant in each processor, and therefore scalable (but not equivalent
+to its serial counterpart). 
+"""
+function FESpaces.ConstantFESpace(
+  model::DistributedDiscreteModel;
+  constraint_type=:global,kwargs...
+)
+  @assert constraint_type ∈ [:global,:local]
+  if constraint_type == :global
+    msg = "ConstantFESpace is NOT scalable in parallel. For testing purposes only."
+    @warn msg
+  end
+
+  spaces = map(local_views(model)) do model
+    ConstantFESpace(model;kwargs...)
+  end
+
+  # Single dof, owned by processor 1 (ghost for all other processors)
+  nranks = length(spaces)
+  cell_gids = get_cell_gids(model)
+  indices = map(partition(cell_gids)) do cell_indices
+    me = part_id(cell_indices)
+    if constraint_type == :global
+      LocalIndices(1,me,Int[1],Int32[1])
+    else
+      LocalIndices(nranks,me,Int[me],Int32[me])
+    end
+  end
+  gids = PRange(indices)
+
+  trian = DistributedTriangulation(map(get_triangulation,spaces),model)
+  vector_type = _find_vector_type(spaces,gids)
+  return DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
+end

src/Geometry.jl

@@ -475,7 +475,8 @@ function Geometry.get_background_model(a::DistributedTriangulation)
 end
 
 function Geometry.num_cells(a::DistributedTriangulation{Df}) where Df
-  gids = get_face_gids(a.model,Df)
+  model = get_background_model(a)
+  gids = get_face_gids(model,Df)
   n_loc_ocells = map(local_views(a),partition(gids)) do a, gids
     glue = get_glue(a,Val(Df))
     @assert isa(glue,FaceToFaceGlue)
@@ -504,11 +505,21 @@ function Geometry.BoundaryTriangulation(
   BoundaryTriangulation(no_ghost,model;kwargs...)
 end
 
+function Geometry.BoundaryTriangulation(
+  trian::DistributedTriangulation;kwargs...)
+  BoundaryTriangulation(no_ghost,trian;kwargs...)
+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.Triangulation(
   portion,::Type{ReferenceFE{Dt}},model::DistributedDiscreteModel{Dm};kwargs...) where {Dt,Dm}
   # Generate global ordering for the faces of dimension Dt (if needed)
@@ -537,6 +548,33 @@ function Geometry.SkeletonTriangulation(
   DistributedTriangulation(trians,model)
 end
 
+# NOTE: The following constructors require adding back the ghost cells: 
+# Potentially, the input `trian` has had some/all of its ghost cells removed. If we do not
+# add them back, some skeleton facets might look like boundary facets to the local constructors...
+function Geometry.BoundaryTriangulation(
+  portion,trian::DistributedTriangulation;kwargs...
+)
+  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...)
+  end
+  DistributedTriangulation(trians,model)
+end
+
+function Geometry.SkeletonTriangulation(
+  portion,trian::DistributedTriangulation;kwargs...
+)
+  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...)
+  end
+  DistributedTriangulation(trians,model)
+end
+
 function Geometry.Triangulation(
   portion,gids::AbstractLocalIndices, args...;kwargs...)
   trian = Triangulation(args...;kwargs...)
@@ -562,8 +600,8 @@ function Geometry.InterfaceTriangulation(
 end
 
 function Geometry.InterfaceTriangulation(a::DistributedTriangulation,b::DistributedTriangulation)
-  trians = map(InterfaceTriangulation,a.trians,b.trians)
   @assert a.model === b.model
+  trians = map(InterfaceTriangulation,a.trians,b.trians)
   DistributedTriangulation(trians,a.model)
 end
 
@@ -623,15 +661,19 @@ function remove_ghost_cells(trian::Triangulation,gids)
   remove_ghost_cells(glue,trian,gids)
 end
 
-function remove_ghost_cells(trian::Union{SkeletonTriangulation,BoundaryTriangulation},gids)
+function remove_ghost_cells(
+  trian::Union{SkeletonTriangulation,BoundaryTriangulation,Geometry.CompositeTriangulation},
+  gids
+)
   model = get_background_model(trian)
   Dm    = num_cell_dims(model)
   glue  = get_glue(trian,Val(Dm))
   remove_ghost_cells(glue,trian,gids)
 end
 
 function remove_ghost_cells(
-  trian::AdaptedTriangulation{Dc,Dp,<:Union{SkeletonTriangulation,BoundaryTriangulation}},gids) where {Dc,Dp}
+  trian::AdaptedTriangulation{Dc,Dp,<:Union{SkeletonTriangulation,BoundaryTriangulation}},gids
+) where {Dc,Dp}
   remove_ghost_cells(trian.trian,gids)
 end
 
@@ -679,7 +721,7 @@ function _find_owned_skeleton_facets(glue,gids)
 end
 
 function add_ghost_cells(dtrian::DistributedTriangulation)
-  dmodel = dtrian.model
+  dmodel = get_background_model(dtrian)
   add_ghost_cells(dmodel,dtrian)
 end
 
@@ -729,7 +771,7 @@ function add_ghost_cells(
 end
 
 function generate_cell_gids(dtrian::DistributedTriangulation)
-  dmodel = dtrian.model
+  dmodel = get_background_model(dtrian)
   generate_cell_gids(dmodel,dtrian)
 end
 

test/AdaptivityCartesianTests.jl

@@ -152,7 +152,7 @@ function main(distribute,ncells,isperiodic)
   fine_adaptivity_glue = get_adaptivity_glue(redist_child_1)
 
   # Redistribute by dispatching on the DistributedCartesianDescriptor
-  pdesc = redist_child_1.metadata
+  pdesc = redist_child_1.metadata.model_metadata
   redist_child_2, redist_glue_child = redistribute(child,pdesc)
 
   # Tests

test/ConstantFESpacesTests.jl

@@ -0,0 +1,38 @@
+module ConstantFESpacesTests
+
+using Test
+using Gridap
+using GridapDistributed, PartitionedArrays
+
+using Gridap.FESpaces, Gridap.Arrays, Gridap.Algebra
+
+function main(distribute,np)
+  @assert prod(np) == 4
+  ranks = distribute(LinearIndices((prod(np),)))
+
+  model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(10,10))
+
+  Ω = Triangulation(model)
+  dΩ = Measure(Ω,2)
+
+  reffe = ReferenceFE(lagrangian, Float64, 1)
+  U = FESpace(model,reffe)
+  V = ConstantFESpace(model)
+  X = MultiFieldFESpace([U,V])
+
+  gids = get_free_dof_ids(V)
+  @test length(gids) == 1
+
+  a((u,λ),(v,μ)) = ∫(u*v + λ*v + μ*u)dΩ 
+  A = assemble_matrix(a,X,X)
+
+  V2 = ConstantFESpace(model;constraint_type=:local)
+  X2 = MultiFieldFESpace([U,V2])
+
+  gids = get_free_dof_ids(V2)
+  @test length(gids) == 4
+
+  A2 = assemble_matrix(a,X2,X2)
+end
+
+end

test/GeometryTests.jl

@@ -96,14 +96,19 @@ function main(distribute,parts)
     add_tag!(labels,"fluid",[fluid])
     cell_to_entity
   end
-  cell_gids=get_cell_gids(model)
-  vcache=PartitionedArrays.p_vector_cache(cell_to_entity,partition(cell_gids))
-  assemble!((a,b)->b, cell_to_entity, map(reverse,vcache) ) |> wait # Make tags consistent
+  cell_gids = get_cell_gids(model)
+  consistent!(PVector(cell_to_entity,partition(cell_gids))) |> wait # Make tags consistent
 
   Ωs = Interior(model,tags="solid")
   Ωf = Interior(model,tags="fluid")
   Γfs = Interface(Ωf,Ωs)
 
+  # CompositeTriangulations
+  Γf = Boundary(Ωf)
+  Λf = Skeleton(Ωf)
+  Λs = Skeleton(Ωs)
+  Γs = Boundary(Ωs)
+
 end
 
 function test_local_part_face_labelings_consistency(lmodel::CartesianDiscreteModel{D},gids,gmodel) where {D}

test/TestApp/src/TestApp.jl

@@ -19,4 +19,5 @@ module TestApp
   include("../../BlockSparseMatrixAssemblersTests.jl")
   include("../../VisualizationTests.jl")
   include("../../AutodiffTests.jl")
+  include("../../ConstantFESpacesTests.jl")
 end