From 33c60cac47cac8b8a33d77b8a9383901a113eed2 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 22 Sep 2015 18:15:54 +0200 Subject: [PATCH] Compute multigrid dof handler also for distributed policy 3D --- include/deal.II/dofs/dof_accessor.templates.h | 28 ++++- source/dofs/dof_handler_policy.cc | 108 ++++++++++++++++-- 2 files changed, 124 insertions(+), 12 deletions(-) diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index d405829cb2..f879e4e0f9 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -1433,12 +1433,30 @@ namespace internal template static - void set_mg_dof_indices (const dealii::DoFAccessor<3, DH,level_dof_access> &, - const int, - const std::vector &, - const unsigned int) + void set_mg_dof_indices (const dealii::DoFAccessor<3, DH,level_dof_access> &accessor, + const int level, + const std::vector &dof_indices, + const unsigned int fe_index) { - AssertThrow (false, ExcNotImplemented ()); //TODO[TH]: implement + const FiniteElement &fe = accessor.get_dof_handler ().get_fe ()[fe_index]; + std::vector::const_iterator next = dof_indices.begin (); + + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex) + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + accessor.set_mg_vertex_dof_index(level, vertex, dof, *next++, fe_index); + + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell; ++line) + for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof) + accessor.line(line)->set_mg_dof_index(level, dof, *next++); + + for (unsigned int quad = 0; quad < GeometryInfo<3>::quads_per_cell; ++quad) + for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof) + accessor.quad(quad)->set_mg_dof_index(level, dof, *next++); + + for (unsigned int dof = 0; dof < fe.dofs_per_hex; ++dof) + accessor.set_mg_dof_index(level, dof, *next++); + + Assert (next == dof_indices.end (), ExcInternalError ()); } }; diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc index dd3eaff84e..7dd0650180 100644 --- a/source/dofs/dof_handler_policy.cc +++ b/source/dofs/dof_handler_policy.cc @@ -845,14 +845,108 @@ namespace internal template static void - renumber_mg_dofs (const std::vector &, - const IndexSet &, - DoFHandler<3,spacedim> &, - const unsigned int , - const bool ) + renumber_mg_dofs (const std::vector &new_numbers, + const IndexSet &indices, + DoFHandler<3,spacedim> &dof_handler, + const unsigned int level, + const bool check_validity) { - // TODO - AssertThrow(false, ExcNotImplemented()); + if (level>=dof_handler.get_tria().n_levels()) + return; + for (typename std::vector::MGVertexDoFs>::iterator i=dof_handler.mg_vertex_dofs.begin(); + i!=dof_handler.mg_vertex_dofs.end(); ++i) + // if the present vertex lives on + // the present level + if ((i->get_coarsest_level() <= level) && + (i->get_finest_level() >= level)) + for (unsigned int d=0; dget_index (level, d); + if (idx != DoFHandler<3>::invalid_dof_index) + i->set_index (level, d, + ((indices.n_elements() == 0) ? + new_numbers[idx] : + new_numbers[indices.index_within_set(idx)])); + + if (check_validity) + Assert(idx != DoFHandler<3>::invalid_dof_index, ExcInternalError ()); + } + + if (dof_handler.get_fe().dofs_per_line > 0 || + dof_handler.get_fe().dofs_per_quad > 0) + { + // save user flags as they will be modified + std::vector user_flags; + dof_handler.get_tria().save_user_flags(user_flags); + const_cast &>(dof_handler.get_tria()).clear_user_flags (); + + // flag all lines adjacent to cells of the current level, as + // those lines logically belong to the same level as the cell, + // at least for isotropic refinement + for (typename DoFHandler<3,spacedim>::level_cell_iterator cell = dof_handler.begin(level); + cell != dof_handler.end(level); ++cell) + for (unsigned int line=0; line < GeometryInfo<3>::lines_per_cell; ++line) + cell->line(line)->set_user_flag(); + + for (typename DoFHandler<3,spacedim>::cell_iterator cell = dof_handler.begin(); + cell != dof_handler.end(); ++cell) + for (unsigned int l=0; l::lines_per_cell; ++l) + if (cell->line(l)->user_flag_set()) + { + for (unsigned int d=0; dline(l)->mg_dof_index(level, d); + if (idx != DoFHandler<3>::invalid_dof_index) + cell->line(l)->set_mg_dof_index (level, d, ((indices.n_elements() == 0) ? + new_numbers[idx] : + new_numbers[indices.index_within_set(idx)])); + if (check_validity) + Assert(idx != DoFHandler<3>::invalid_dof_index, ExcInternalError ()); + } + cell->line(l)->clear_user_flag(); + } + + // flag all quads adjacent to cells of the current level, as + // those quads logically belong to the same level as the cell, + // at least for isotropic refinement + for (typename DoFHandler<3,spacedim>::level_cell_iterator cell = dof_handler.begin(level); + cell != dof_handler.end(level); ++cell) + for (unsigned int quad=0; quad < GeometryInfo<3>::quads_per_cell; ++quad) + cell->quad(quad)->set_user_flag(); + + for (typename DoFHandler<3,spacedim>::cell_iterator cell = dof_handler.begin(); + cell != dof_handler.end(); ++cell) + for (unsigned int l=0; l::quads_per_cell; ++l) + if (cell->quad(l)->user_flag_set()) + { + for (unsigned int d=0; dquad(l)->mg_dof_index(level, d); + if (idx != DoFHandler<1>::invalid_dof_index) + cell->quad(l)->set_mg_dof_index (level, d, ((indices.n_elements() == 0) ? + new_numbers[idx] : + new_numbers[indices.index_within_set(idx)])); + if (check_validity) + Assert(idx != DoFHandler<3>::invalid_dof_index, ExcInternalError ()); + } + cell->quad(l)->clear_user_flag(); + } + + // finally, restore user flags + const_cast &>(dof_handler.get_tria()).load_user_flags (user_flags); + } + + for (std::vector::iterator i=dof_handler.mg_levels[level]->dof_object.dofs.begin(); + i!=dof_handler.mg_levels[level]->dof_object.dofs.end(); ++i) + { + if (*i != DoFHandler<3>::invalid_dof_index) + { + Assert(*i