From 0ea6c67e8aa54094b799c2b6c0739424324bd125 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 8 Nov 2017 17:24:00 -0700 Subject: [PATCH] Generalize the numbering in the Cuthill-McKee algorithm. --- include/deal.II/dofs/dof_renumbering.h | 67 +++++++++---- source/dofs/dof_renumbering.cc | 133 +++++++++++++++++++++---- 2 files changed, 161 insertions(+), 39 deletions(-) diff --git a/include/deal.II/dofs/dof_renumbering.h b/include/deal.II/dofs/dof_renumbering.h index c61617cbc4..99ce70410a 100644 --- a/include/deal.II/dofs/dof_renumbering.h +++ b/include/deal.II/dofs/dof_renumbering.h @@ -518,15 +518,6 @@ namespace DoFRenumbering * comparison of various algorithms in the documentation of the * DoFRenumbering namespace. * - * If the given DoFHandler uses a distributed triangulation (i.e., if - * dof_handler.locally_owned() is not the complete index set), the - * renumbering is performed on each processor's degrees of freedom - * individually, without any communication between processors. In other - * words, the resulting renumbering is an attempt at minimizing the bandwidth - * of each diagonal block of the matrix corresponding to one processor - * separately, without making an attempt at minimizing the bandwidth of - * the global matrix. - * * @param dof_handler The DoFHandler or hp::DoFHandler object to work on. * @param reversed_numbering Whether to use the original Cuthill-McKee * algorithm, or to reverse the ordering. @@ -535,14 +526,50 @@ namespace DoFRenumbering * @param starting_indices A set of degrees of freedom that form the first * level of renumbered degrees of freedom. If the set is empty, then a * single starting entry is chosen automatically among those that have the - * smallest number of others that couple with it. If the DoFHandler is built - * on a parallel triangulation, then on every processor, these starting - * indices need to be a (possibly empty) subset of the - * @ref GlossLocallyOwnedDof "locally owned degrees of freedom". - * These will then be used as starting indices for the local renumbering on - * the current processor. (In other words, you will have to choose this - * argument differently on every processor, unless of course you pass an - * empty list as is the default.) + * smallest number of others that couple with it. + * + *

Operation in parallel

+ * + * If the given DoFHandler uses a distributed triangulation (i.e., if + * dof_handler.locally_owned() is not the complete index set), the + * renumbering is performed on each processor's degrees of freedom + * individually, without any communication between processors. In other + * words, the resulting renumbering is an attempt at minimizing the bandwidth + * of each diagonal block of the matrix corresponding to one processor + * separately, without making an attempt at minimizing the bandwidth of + * the global matrix. Furthermore, the renumbering reuses exactly the + * same set of DoF indices that each processor used before. In other words, + * if the previous numbering of DoFs on one processor used a contiguous + * range of DoF indices, then so will the DoFs on that processor after + * the renumbering, and they will occupy the same range. The same is true + * if the previous numbering of DoFs on a processor consisted of a number + * of index ranges or single indices: after renumbering, the locally owned + * DoFs on that processor will use the exact same indices, just in a + * different order. + * + * In addition, if the DoFHandler is built on a parallel triangulation, then + * on every processor, the starting indices for renumbering need to be a + * (possibly empty) subset of the + * @ref GlossLocallyActiveDof "locally active degrees of freedom". In + * general, these starting indices will be different on each processor + * (unless of course you pass an empty list as is the default), + * and each processor will use them as starting indices for the local + * renumbering on that processor. + * + * The starting indices must be locally active degrees of freedom, but + * the function will only renumber the locally owned subset of the + * locally owned DoFs. The function accepts starting indices from the + * largest set of locally active degrees of freedom because a typical + * renumbering operation with this function starts with indices that + * are located on the boundary -- in the case of the current function, + * that would be the boundary between processor subdomains. Since the + * degrees of freedom that are located on subdomain interfaces may + * be owned by either one of the two processors that own the adjacent + * subdomains, it is not always easy to identify starting indices that + * are locally owned. On the other hand, all degrees of freedom on subdomain + * interfaces are locally active, and so the function accepts them as + * starting indices even though it can only renumber them on a given + * processor if they are also locally owned. */ template void @@ -554,8 +581,10 @@ namespace DoFRenumbering /** * Compute the renumbering vector needed by the Cuthill_McKee() function. - * Does not perform the renumbering on the DoFHandler dofs but returns the - * renumbering vector. + * This function does not perform the renumbering on the DoFHandler DoFs but + * only returns the renumbering vector. + * + * See the Cuthill_McKee() function for an explanation of the arguments. */ template void diff --git a/source/dofs/dof_renumbering.cc b/source/dofs/dof_renumbering.cc index dbdd1364da..df951276ba 100644 --- a/source/dofs/dof_renumbering.cc +++ b/source/dofs/dof_renumbering.cc @@ -417,27 +417,55 @@ namespace DoFRenumbering // local index space, i.e., the locally owned part of the // sparsity pattern. // - // create first the global sparsity pattern, and then the local + // first figure out whether the user only gave us starting + // indices that are locally owned, or that are only locally + // relevant. in the process, also check that all indices + // really belong to at least the locally relevant ones + IndexSet locally_active_dofs; + DoFTools::extract_locally_active_dofs(dof_handler, locally_active_dofs); + + bool needs_locally_active = false; + for (unsigned int i=0; i row_entries; - for (unsigned int i=0; i local_starting_indices (starting_indices.size()); for (unsigned int i=0; i my_new_indices (index_set_to_use.n_elements()); + SparsityTools::reorder_Cuthill_McKee (local_sparsity, my_new_indices, local_starting_indices); if (reversed_numbering) - new_indices = Utilities::reverse_permutation (new_indices); + my_new_indices = Utilities::reverse_permutation (my_new_indices); + + // now that we have a re-enumeration of all DoFs, we need to throw + // out the ones that are not locally owned in case we have worked + // with the locally active ones. that's because the renumbering + // functions only want new indices for the locally owned DoFs (other + // processors are responsible for renumbering the ones that are + // on cell interfaces) + if (needs_locally_active == true) + { + // first step: figure out which DoF indices to eliminate + IndexSet active_but_not_owned_dofs = locally_active_dofs; + active_but_not_owned_dofs.subtract_set(locally_owned_dofs); + + std::set erase_these_indices; + for (const auto &p : active_but_not_owned_dofs) + { + const auto index = index_set_to_use.index_within_set(p); + Assert (index < index_set_to_use.n_elements(), ExcInternalError()); + erase_these_indices.insert (my_new_indices[index]); + my_new_indices[index] = numbers::invalid_dof_index; + } + Assert (erase_these_indices.size() == active_but_not_owned_dofs.n_elements(), + ExcInternalError()); + Assert (std::count (my_new_indices.begin(), + my_new_indices.end(), + numbers::invalid_dof_index) == + active_but_not_owned_dofs.n_elements(), + ExcInternalError()); + + // then compute a renumbering of the remaining ones + std::vector translate_indices (my_new_indices.size()); + { + std::set::const_iterator + next_erased_index = erase_these_indices.begin(); + types::global_dof_index next_new_index = 0; + for (unsigned int i=0; i