GMG working in serial

src/LinearSolvers/GMGLinearSolvers.jl

@@ -415,7 +415,9 @@ function gmg_work_vectors(smatrices::AbstractVector{<:AbstractMatrix})
   return work_vectors
 end
 
-function apply_GMG_level!(lev::Integer,xh::Union{PVector,Nothing},rh::Union{PVector,Nothing},ns::GMGNumericalSetup)
+function apply_GMG_level!(
+  lev::Integer,xh::Union{<:AbstractVector,Nothing},rh::Union{<:AbstractVector,Nothing},ns::GMGNumericalSetup
+)
   mh = ns.solver.mh
   parts = get_level_parts(mh,lev)
   if i_am_in(parts)

src/MultilevelTools/DistributedGridTransferOperators.jl

@@ -65,11 +65,11 @@ function _get_interpolation_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mod
   if i_am_in(cparts)
     model_h = get_model_before_redist(sh,lev)
     Uh   = get_fe_space_before_redist(sh,lev)
-    fv_h = pfill(0.0,partition(Uh.gids))
+    fv_h = zero_free_values(Uh)
     dv_h = (mode == :solution) ? get_dirichlet_dof_values(Uh) : zero_dirichlet_values(Uh)
 
     UH   = get_fe_space(sh,lev+1)
-    fv_H = pfill(0.0,partition(UH.gids))
+    fv_H = zero_free_values(UH)
     dv_H = (mode == :solution) ? get_dirichlet_dof_values(UH) : zero_dirichlet_values(UH)
 
     cache_refine = model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H
@@ -222,7 +222,7 @@ end
 ### Applying the operators:
 
 # A) Prolongation, without redistribution
-function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:prolongation},Val{false}},x::PVector)
+function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:prolongation},Val{false}},x::AbstractVector)
   cache_refine, cache_redist = A.cache
   model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H = cache_refine
 
@@ -235,7 +235,7 @@ function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:p
 end
 
 # B.1) Restriction, without redistribution, by interpolation
-function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:interpolation}},x::PVector)
+function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:interpolation}},x::AbstractVector)
   cache_refine, cache_redist = A.cache
   model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H = cache_refine
 
@@ -248,7 +248,7 @@ function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:r
 end
 
 # B.2) Restriction, without redistribution, by projection
-function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:projection}},x::PVector)
+function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:projection}},x::AbstractVector)
   cache_refine, cache_redist = A.cache
   model_h, Uh, fv_h, dv_h, VH, AH_ns, lH, xH, bH, bH0, assem = cache_refine
 
@@ -265,7 +265,7 @@ function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:r
 end
 
 # B.3) Restriction, without redistribution, by dof selection (only nodal dofs)
-function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:dof_mask}},x::PVector)
+function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:dof_mask}},x::AbstractVector)
   cache_refine, cache_redist = A.cache
   model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H = cache_refine
 
@@ -367,6 +367,19 @@ end
 
 ###############################################################
 
+function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:dual_projection}},x::AbstractVector)
+  cache_refine, cache_redist = A.cache
+  model_h, Uh, VH, Mh_ns, rh, uh, assem, dΩhH = cache_refine
+  fv_h = get_free_dof_values(uh)
+
+  solve!(rh,Mh_ns,x)
+  copy!(fv_h,rh)
+  v = get_fe_basis(VH)
+  assemble_vector!(y,assem,collect_cell_vector(VH,∫(v⋅uh)*dΩhH))
+
+  return y
+end
+
 function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:dual_projection}},x::PVector)
   cache_refine, cache_redist = A.cache
   model_h, Uh, VH, Mh_ns, rh, uh, assem, dΩhH = cache_refine

src/MultilevelTools/FESpaceHierarchies.jl

@@ -56,6 +56,16 @@ function _cell_conformity(
   return CellConformity(cell_reffe,conformity)
 end
 
+function _cell_conformity(
+  model::DiscreteModel,
+  reffe::ReferenceFE; 
+  conformity=nothing, kwargs...
+) :: CellConformity
+  cell_reffe = Fill(reffe,num_cells(model))
+  conformity = Conformity(Gridap.Arrays.testitem(cell_reffe),conformity)
+  return CellConformity(cell_reffe,conformity)
+end
+
 function _cell_conformity(
   model::GridapDistributed.DistributedDiscreteModel,args...;kwargs...
 ) :: AbstractVector{<:CellConformity}

src/MultilevelTools/ModelHierarchies.jl

@@ -139,6 +139,44 @@ function CartesianModelHierarchy(
   return HierarchicalArray(meshes,level_parts)
 end
 
+function ModelHierarchy(models::Vector{<:DiscreteModel})
+  nlevs = length(models)
+  @check all(map(i -> isa(models[i],AdaptedDiscreteModel),1:nlevs-1)) "Hierarchy models are not adapted models."
+  for lev in 1:nlevs-1
+    @check Adaptivity.is_child(models[lev],models[lev+1]) "Incorrect hierarchy of models."
+  end
+
+  level_parts = fill(DebugArray([1]),nlevs)
+  meshes = Vector{ModelHierarchyLevel}(undef,nlevs)
+  for lev in 1:nlevs-1
+    glue  = Gridap.Adaptivity.get_adaptivity_glue(models[lev])
+    model = models[lev]
+    meshes[lev] = ModelHierarchyLevel(lev,model,glue,nothing,nothing)
+  end
+  meshes[nlevs] = ModelHierarchyLevel(nlevs,models[nlevs],nothing,nothing,nothing)
+  return HierarchicalArray(meshes,level_parts)
+end
+
+function ModelHierarchy(models::Vector{<:GridapDistributed.DistributedDiscreteModel})
+  nlevs = length(models)
+  @check all(map(i -> isa(models[i],GridapDistributed.DistributedAdaptedDiscreteModel),1:nlevs-1)) "Hierarchy models are not adapted models."
+  for lev in 1:nlevs-1
+    @check Adaptivity.is_child(models[lev],models[lev+1]) "Incorrect hierarchy of models."
+  end
+  ranks = get_parts(models[1])
+  @check all(m -> get_parts(m) === ranks, models) "Models have different communicators."
+
+  level_parts = fill(ranks,nlevs)
+  meshes = Vector{ModelHierarchyLevel}(undef,nlevs)
+  for lev in 1:nlevs-1
+    glue  = Gridap.Adaptivity.get_adaptivity_glue(models[lev])
+    model = models[lev]
+    meshes[lev] = ModelHierarchyLevel(lev,model,glue,nothing,nothing)
+  end
+  meshes[nlevs] = ModelHierarchyLevel(nlevs,models[nlevs],nothing,nothing,nothing)
+  return HierarchicalArray(meshes,level_parts)
+end
+
 """
     P4estCartesianModelHierarchy(
       ranks,np_per_level,domain,nc::NTuple{D,<:Integer};

src/PatchBasedSmoothers/PatchBasedSmoothers.jl

@@ -22,6 +22,7 @@ export setup_patch_prolongation_operators, setup_patch_restriction_operators
 include("seq/PatchDecompositions.jl")
 include("mpi/PatchDecompositions.jl")
 include("seq/PatchTriangulations.jl")
+include("seq/CoarsePatchDecompositions.jl")
 
 # FESpaces
 include("seq/PatchFESpaces.jl")

src/PatchBasedSmoothers/seq/CoarsePatchDecompositions.jl

@@ -0,0 +1,63 @@
+
+# PatchDecomposition where each patch is given by a coarse cell on the parent mesh
+function CoarsePatchDecomposition(
+  model::Gridap.Adaptivity.AdaptedDiscreteModel{Dc,Dp},
+  patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude(),
+  boundary_tag_names::AbstractArray{String}=["boundary"]
+) where {Dc,Dp}
+  glue = Gridap.Adaptivity.get_adaptivity_glue(model)
+
+  patch_cells = glue.o2n_faces_map
+  patch_cells_overlapped = compute_patch_overlapped_cells(patch_cells)
+  patch_facets = get_coarse_patch_facets(model, patch_cells)
+  patch_cells_faces_on_boundary = compute_patch_cells_faces_on_boundary(
+    model, patch_cells, patch_cells_overlapped,
+    patch_facets, patch_boundary_style, boundary_tag_names
+  )
+
+  return PatchDecomposition{Dc,Dc,Dp}(
+    model, patch_cells, patch_cells_overlapped, patch_cells_faces_on_boundary
+  )
+end
+
+function get_coarse_patch_facets(
+  model::Gridap.Adaptivity.AdaptedDiscreteModel{Dc,Dp},patch_cells
+) where {Dc,Dp}
+  topo = get_grid_topology(model)
+  c2f_map = Geometry.get_faces(topo,Dc,Dc-1)
+  f2c_map = Geometry.get_faces(topo,Dc-1,Dc)
+
+  touched = fill(false,num_faces(topo,Dc-1))
+  ptrs = fill(0,length(patch_cells)+1)
+  for (patch,cells) in enumerate(patch_cells)
+    for c in cells
+      for f in c2f_map[c]
+        nbors = f2c_map[f]
+        if !touched[f] && (length(nbors) == 2) && all(n -> n ∈ cells, nbors)
+          touched[f] = true
+          ptrs[patch+1] += 1
+        end
+      end
+    end
+  end
+  Arrays.length_to_ptrs!(ptrs)
+
+  data = fill(0,ptrs[end]-1)
+  fill!(touched,false)
+  for (patch,cells) in enumerate(patch_cells)
+    for c in cells
+      for f in c2f_map[c]
+        nbors = f2c_map[f]
+        if !touched[f] && (length(nbors) == 2) && all(n -> n ∈ cells, nbors)
+          touched[f] = true
+          data[ptrs[patch]] = f
+          ptrs[patch] += 1
+        end
+      end
+    end
+  end
+  Arrays.rewind_ptrs!(ptrs)
+
+  patch_facets = Table(data,ptrs)
+  return patch_facets
+end