From 6b81785317d9b7006943b55a1f1392178a6cb1e2 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 20 Sep 2012 13:37:48 +0000 Subject: [PATCH] Add DoFRenumbering::block_wise functions. git-svn-id: https://svn.dealii.org/trunk@26565 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/deal.II/dofs/dof_renumbering.h | 118 ++++++- deal.II/source/dofs/dof_renumbering.cc | 311 +++++++++++++++++- deal.II/source/dofs/dof_renumbering.inst.in | 17 + 3 files changed, 438 insertions(+), 8 deletions(-) diff --git a/deal.II/include/deal.II/dofs/dof_renumbering.h b/deal.II/include/deal.II/dofs/dof_renumbering.h index bb85d16439..6b7d838b1f 100644 --- a/deal.II/include/deal.II/dofs/dof_renumbering.h +++ b/deal.II/include/deal.II/dofs/dof_renumbering.h @@ -144,18 +144,21 @@ DEAL_II_NAMESPACE_OPEN * may be difficult, however, and in many cases will not justify the effort. * * - *

Component-wise numbering

+ *

Component-wise and block-wise numberings

* * For finite elements composed of several base elements using the FESystem * class, or for elements which provide several components themselves, it * may be of interest to sort the DoF indices by component. This will then * bring out the block matrix structure, since otherwise the degrees of freedom * are numbered cell-wise without taking into account that they may belong to - * different components. + * different components. For example, one may want to sort degree of freedom for + * a Stokes discretization so that we first get all velocities and then all + * the pressures so that the resulting matrix naturally decomposes into a + * $2\times 2$ system. * * This kind of numbering may be obtained by calling the * component_wise() function of this class. Since it does not touch - * the order of indices within each, it may be worthwhile to first + * the order of indices within each component, it may be worthwhile to first * renumber using the Cuthill-McKee or a similar algorithm and * afterwards renumbering component-wise. This will bring out the * matrix structure and additionally have a good numbering within each @@ -163,8 +166,26 @@ DEAL_II_NAMESPACE_OPEN * * The component_wise() function allows not only to honor enumeration based on * vector components, but also allows to group together vector components into - * "blocks". See @ref GlossComponent vs @ref GlossBlock for a description of - * the difference. + * "blocks" using a defaulted argument to the various DoFRenumber::component_wise() + * functinos (see @ref GlossComponent vs @ref GlossBlock for a description of + * the difference). The blocks designated through this argument may, but do not + * have to be, equal to the blocks that the finite element reports. For example, + * a typical Stokes element would be + * @code + * FESystem stokes_fe (FE_Q(2), dim, // dim velocities + * FE_Q(1), 1); // one pressure + * @endcode + * This element has dim+1 vector components and equally many + * blocks. However, one may want to consider the velocities as one logical + * block so that all velocity degrees of freedom are enumerated the same + * way, independent of whether they are $x$- or $y$-velocities. This is done, + * for example, in step-20 and step-22 as well as several other tutorial programs. + * + * On the other hand, if you really want to use block structure reported + * by the finite element itself (a case that is often the case if you have + * finite elements that have multiple vector components, e.g. the FE_RaviartThomas + * or FE_Nedelec elements) then you can use the DoFRenumber::block_wise instead + * of the DoFRenumbering::component_wise functions. * * *

Cell-wise numbering

@@ -700,7 +721,7 @@ namespace DoFRenumbering * target component for dofs with * component @p i in the * FESystem. Naming the same - * component more than once is + * target component more than once is * possible and results in a * blocking of several components * into one. This is discussed in @@ -760,7 +781,7 @@ namespace DoFRenumbering template void component_wise (MGDoFHandler& dof_handler, - unsigned int level, + const unsigned int level, const std::vector& target_component = std::vector()); @@ -799,6 +820,89 @@ namespace DoFRenumbering * @} */ + /** + * @name Block-wise numberings + * @{ + */ + + /** + * Sort the degrees of freedom by + * vector block. The + * numbering within each + * block is not touched, so a + * degree of freedom with index + * $i$, belonging to some + * block, and another degree + * of freedom with index $j$ + * belonging to the same + * block will be assigned new + * indices $n(i)$ and $n(j)$ with + * $n(i)n(j)$ if $i>j$. + */ + template + void + block_wise (DoFHandler &dof_handler); + + + /** + * Sort the degrees of freedom by + * block. It does the same + * thing as the above function. + */ + template + void + block_wise (hp::DoFHandler &dof_handler); + + /** + * Sort the degrees of freedom by + * block. It does the same + * thing as the above function, + * only that it does this for one + * single level of a multi-level + * discretization. The + * non-multigrid part of the + * MGDoFHandler is not touched. + */ + template + void + block_wise (MGDoFHandler &dof_handler, + const unsigned int level); + + + /** + * Sort the degrees of freedom by + * block. It does the same + * thing as the previous + * functions, but more: it + * renumbers not only every level + * of the multigrid part, but + * also the global, + * i.e. non-multigrid components. + */ + template + void + block_wise (MGDoFHandler &dof_handler); + + /** + * Computes the renumbering + * vector needed by the + * block_wise() + * functions. Does not perform + * the renumbering on the + * DoFHandler dofs but returns + * the renumbering vector. + */ + template + unsigned int + compute_block_wise (std::vector& new_dof_indices, + const ITERATOR& start, + const ENDITERATOR& end); + + /** + * @} + */ + /** * @name Various cell-wise numberings * @{ diff --git a/deal.II/source/dofs/dof_renumbering.cc b/deal.II/source/dofs/dof_renumbering.cc index 99bec3b823..01eaf7e9f3 100644 --- a/deal.II/source/dofs/dof_renumbering.cc +++ b/deal.II/source/dofs/dof_renumbering.cc @@ -781,7 +781,7 @@ namespace DoFRenumbering } // now we've got all indices sorted - // into buckets labelled with their + // into buckets labeled by their // target component number. we've // only got to traverse this list // and assign the new indices @@ -899,6 +899,315 @@ namespace DoFRenumbering return next_free_index; } + + + template + void + block_wise (DoFHandler &dof_handler) + { + std::vector renumbering (dof_handler.n_locally_owned_dofs(), + DoFHandler::invalid_dof_index); + + typedef + internal::WrapDoFIterator + ::active_cell_iterator> + ITERATOR; + + typename DoFHandler::active_cell_iterator + istart = dof_handler.begin_active(); + ITERATOR start = istart; + const typename DoFHandler::cell_iterator + end = dof_handler.end(); + + const unsigned int result = + compute_block_wise::cell_iterator> + (renumbering, start, end); + if (result == 0) + return; + + // verify that the last numbered + // degree of freedom is either + // equal to the number of degrees + // of freedom in total (the + // sequential case) or in the + // distributed case at least + // makes sense + Assert ((result == dof_handler.n_locally_owned_dofs()) + || + ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) + && + (result <= dof_handler.n_dofs())), + ExcRenumberingIncomplete()); + + dof_handler.renumber_dofs (renumbering); + } + + + + template + void + block_wise (hp::DoFHandler &dof_handler) + { + std::vector renumbering (dof_handler.n_dofs(), + hp::DoFHandler::invalid_dof_index); + + typedef + internal::WrapDoFIterator::active_cell_iterator> ITERATOR; + + typename hp::DoFHandler::active_cell_iterator + istart = dof_handler.begin_active(); + ITERATOR start = istart; + const typename hp::DoFHandler::cell_iterator + end = dof_handler.end(); + + const unsigned int result = + compute_block_wise::cell_iterator>(renumbering, + start, end); + + if (result == 0) + return; + + Assert (result == dof_handler.n_dofs(), + ExcRenumberingIncomplete()); + + dof_handler.renumber_dofs (renumbering); + } + + + + template + void + block_wise (MGDoFHandler &dof_handler, + const unsigned int level) + { + std::vector renumbering (dof_handler.n_dofs(level), + DoFHandler::invalid_dof_index); + + typedef + internal::WrapMGDoFIterator::cell_iterator> ITERATOR; + + typename MGDoFHandler::cell_iterator + istart =dof_handler.begin(level); + ITERATOR start = istart; + typename MGDoFHandler::cell_iterator + iend = dof_handler.end(level); + const ITERATOR end = iend; + + const unsigned int result = + compute_block_wise( + renumbering, start, end); + + if (result == 0) return; + + Assert (result == dof_handler.n_dofs(level), + ExcRenumberingIncomplete()); + + if (renumbering.size()!=0) + dof_handler.renumber_dofs (level, renumbering); + } + + + + template + void + block_wise (MGDoFHandler &dof_handler) + { + // renumber the non-MG part of + // the DoFHandler in parallel to + // the MG part. Because + // MGDoFHandler::renumber_dofs + // uses the user flags we can't + // run renumbering on individual + // levels in parallel to the + // other levels + void (*non_mg_part) (DoFHandler &) + = &block_wise; + Threads::Task<> + task = Threads::new_task (non_mg_part, dof_handler); + + for (unsigned int level=0; level + unsigned int + compute_block_wise (std::vector& new_indices, + const ITERATOR & start, + const ENDITERATOR& end) + { + const hp::FECollection + fe_collection (start->get_dof_handler().get_fe ()); + + // do nothing if the FE has only + // one component + if (fe_collection.n_blocks() == 1) + { + new_indices.resize(0); + return 0; + } + + // vector to hold the dof indices on + // the cell we visit at a time + std::vector local_dof_indices; + + // prebuilt list to which block + // a given dof on a cell + // should go. + std::vector > block_list (fe_collection.size()); + for (unsigned int f=0; f & fe = fe_collection[f]; + block_list[f].resize(fe.dofs_per_cell); + for (unsigned int i=0; i > + block_to_dof_map (fe_collection.n_blocks()); + for (ITERATOR cell=start; cell!=end; ++cell) + if (cell->is_locally_owned()) + { + // on each cell: get dof indices + // and insert them into the global + // list using their component + const unsigned int fe_index = cell->active_fe_index(); + const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell; + local_dof_indices.resize (dofs_per_cell); + cell.get_dof_indices (local_dof_indices); + for (unsigned int i=0; iget_dof_handler().locally_owned_dofs().is_element(local_dof_indices[i])) + block_to_dof_map[block_list[fe_index][i]]. + push_back (local_dof_indices[i]); + } + + // now we've got all indices sorted + // into buckets labeled by their + // target block number. we've + // only got to traverse this list + // and assign the new indices + // + // however, we first want to sort + // the indices entered into the + // buckets to preserve the order + // within each component and during + // this also remove duplicate + // entries + for (unsigned int block=0; block shifts(n_buckets); + + if (const parallel::distributed::Triangulation * tria + = (dynamic_cast*> + (&start->get_dof_handler().get_tria()))) + { +#ifdef DEAL_II_USE_P4EST + std::vector local_dof_count(n_buckets); + + for (unsigned int c=0; c + all_dof_counts(fe_collection.n_components() * + Utilities::MPI::n_mpi_processes (tria->get_communicator())); + + MPI_Allgather ( &local_dof_count[0], n_buckets, MPI_UNSIGNED, &all_dof_counts[0], + n_buckets, MPI_UNSIGNED, tria->get_communicator()); + + for (unsigned int i=0; ilocally_owned_subdomain()+i] + == + local_dof_count[i], + ExcInternalError()); + + //calculate shifts + unsigned int cumulated = 0; + for (unsigned int c=0; clocally_owned_subdomain(); ++i) + shifts[c] += all_dof_counts[c+n_buckets*i]; + for (unsigned int i=0; iget_communicator()); ++i) + cumulated += all_dof_counts[c+n_buckets*i]; + } +#else + (void)tria; + Assert (false, ExcInternalError()); +#endif + } + else + { + shifts[0] = 0; + for (unsigned int c=1; c::const_iterator + begin_of_component = block_to_dof_map[block].begin(), + end_of_component = block_to_dof_map[block].end(); + + next_free_index = shifts[block]; + + for (typename std::vector::const_iterator + dof_index = begin_of_component; + dof_index != end_of_component; ++dof_index) + { + Assert (start->get_dof_handler().locally_owned_dofs() + .index_within_set(*dof_index) + < + new_indices.size(), + ExcInternalError()); + new_indices[start->get_dof_handler().locally_owned_dofs() + .index_within_set(*dof_index)] + = next_free_index++; + } + } + + return next_free_index; + } + + + namespace { // helper function for hierarchical() diff --git a/deal.II/source/dofs/dof_renumbering.inst.in b/deal.II/source/dofs/dof_renumbering.inst.in index 769af6b944..fd9feb3cd4 100644 --- a/deal.II/source/dofs/dof_renumbering.inst.in +++ b/deal.II/source/dofs/dof_renumbering.inst.in @@ -44,6 +44,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS (DoFHandler&, const std::vector&); + template + void block_wise + (DoFHandler&); + template void subdomain_wise > (DoFHandler &); @@ -169,6 +173,19 @@ namespace DoFRenumbering (MGDoFHandler&, const std::vector&); + template + void block_wise + (hp::DoFHandler&); + + template + void block_wise + (MGDoFHandler&, + unsigned int); + + template + void block_wise + (MGDoFHandler&); + template void hierarchical (DoFHandler&); -- 2.39.5