]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Generalize the numbering in the Cuthill-McKee algorithm.
authorWolfgang Bangerth <bangerth@colostate.edu>
Thu, 9 Nov 2017 00:24:00 +0000 (17:24 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 27 Nov 2017 17:40:36 +0000 (10:40 -0700)
include/deal.II/dofs/dof_renumbering.h
source/dofs/dof_renumbering.cc

index c61617cbc459d7d232f0afbd2ba836df9e183e4f..99ce70410a5769fd3c6fdf9ed0841b50ecfeb138 100644 (file)
@@ -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 <i>each diagonal block of the matrix corresponding to one processor</i>
-   * 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.
+   *
+   * <h4> Operation in parallel </h4>
+   *
+   * 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 <i>each diagonal block of the matrix corresponding to one processor</i>
+   * 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 <typename DoFHandlerType>
   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 <typename DoFHandlerType>
   void
index dbdd1364dabd2121dfd80fa3bf0446814f5a617f..df951276ba249a1ad60338ec05ca03a761c818a3 100644 (file)
@@ -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<starting_indices.size(); ++i)
+          {
+            if ((needs_locally_active == /* previously already set to */true)
+                ||
+                (locally_owned_dofs.is_element (starting_indices[i]) == false))
+              {
+                Assert (locally_active_dofs.is_element (starting_indices[i]),
+                        ExcMessage ("You specified global degree of freedom "
+                                    + Utilities::to_string(starting_indices[i]) +
+                                    " as a starting index, but this index is not among the "
+                                    "locally active ones on this processor, as required "
+                                    "for this function."));
+                needs_locally_active = true;
+              }
+          }
+
+        const IndexSet index_set_to_use = (needs_locally_active ?
+                                           locally_active_dofs :
+                                           locally_owned_dofs);
+
+        // then create first the global sparsity pattern, and then the local
         // sparsity pattern from the global one by transferring its indices to
-        // processor-local index space
+        // processor-local (locally owned or locally active) index space
         DynamicSparsityPattern dsp (dof_handler.n_dofs(),
                                     dof_handler.n_dofs(),
-                                    locally_owned_dofs);
+                                    index_set_to_use);
         DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints);
 
-        DynamicSparsityPattern local_sparsity(locally_owned_dofs.n_elements(),
-                                              locally_owned_dofs.n_elements());
+        DynamicSparsityPattern local_sparsity(index_set_to_use.n_elements(),
+                                              index_set_to_use.n_elements());
         std::vector<types::global_dof_index> row_entries;
-        for (unsigned int i=0; i<locally_owned_dofs.n_elements(); ++i)
+        for (unsigned int i=0; i<index_set_to_use.n_elements(); ++i)
           {
-            const types::global_dof_index row = locally_owned_dofs.nth_index_in_set(i);
+            const types::global_dof_index row = index_set_to_use.nth_index_in_set(i);
             const unsigned int row_length = dsp.row_length(row);
             row_entries.clear();
             for (unsigned int j=0; j<row_length; ++j)
               {
                 const unsigned int col = dsp.column_number(row, j);
-                if (col != row && locally_owned_dofs.is_element(col))
-                  row_entries.push_back(locally_owned_dofs.index_within_set(col));
+                if (col != row && index_set_to_use.is_element(col))
+                  row_entries.push_back(index_set_to_use.index_within_set(col));
               }
             local_sparsity.add_entries(i, row_entries.begin(), row_entries.end(),
                                        true);
@@ -446,23 +474,88 @@ namespace DoFRenumbering
         // translate starting indices from global to local indices
         std::vector<types::global_dof_index> local_starting_indices (starting_indices.size());
         for (unsigned int i=0; i<starting_indices.size(); ++i)
-          {
-            Assert (locally_owned_dofs.is_element (starting_indices[i]),
-                    ExcMessage ("You specified global degree of freedom "
-                                + Utilities::to_string(starting_indices[i]) +
-                                " as a starting index, but this index is not among the "
-                                "locally owned ones on this processor."));
-            local_starting_indices[i] = locally_owned_dofs.index_within_set(starting_indices[i]);
-          }
+          local_starting_indices[i] = index_set_to_use.index_within_set(starting_indices[i]);
 
         // then do the renumbering on the locally owned portion
         AssertDimension(new_indices.size(), locally_owned_dofs.n_elements());
-        SparsityTools::reorder_Cuthill_McKee (local_sparsity, new_indices,
+        std::vector<types::global_dof_index> 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<types::global_dof_index> 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<types::global_dof_index> translate_indices (my_new_indices.size());
+            {
+              std::set<types::global_dof_index>::const_iterator
+              next_erased_index = erase_these_indices.begin();
+              types::global_dof_index next_new_index = 0;
+              for (unsigned int i=0; i<translate_indices.size(); ++i)
+                if ((next_erased_index != erase_these_indices.end())
+                    &&
+                    (*next_erased_index == i))
+                  {
+                    translate_indices[i] = numbers::invalid_dof_index;
+                    ++next_erased_index;
+                  }
+                else
+                  {
+                    translate_indices[i] = next_new_index;
+                    ++next_new_index;
+                  }
+              Assert (next_new_index == locally_owned_dofs.n_elements(),
+                      ExcInternalError());
+            }
+
+            // and then do the renumbering of the result of the
+            // Cuthill-McKee algorithm above, right into the output array
+            new_indices.clear();
+            new_indices.reserve(locally_owned_dofs.n_elements());
+            for (const auto &p : my_new_indices)
+              if (p != numbers::invalid_dof_index)
+                {
+                  Assert (translate_indices[p] != numbers::invalid_dof_index,
+                          ExcInternalError());
+                  new_indices.push_back (translate_indices[p]);
+                }
+            Assert (new_indices.size() == locally_owned_dofs.n_elements(),
+                    ExcInternalError());
+          }
+        else
+          new_indices = std::move (my_new_indices);
 
-        // convert indices back to global index space
+        // convert indices back to global index space. in both of the branches
+        // above, we ended up with new_indices only containing the local
+        // indices of the locally-owned DoFs. so that's where we get the
+        // indices
         for (std::size_t i=0; i<new_indices.size(); ++i)
           new_indices[i] = locally_owned_dofs.nth_index_in_set(new_indices[i]);
       }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.