]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Split out the function that does the actual Cuthill-McKee algorithm from DoFRenumberi...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 31 Oct 2008 17:24:43 +0000 (17:24 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 31 Oct 2008 17:24:43 +0000 (17:24 +0000)
git-svn-id: https://svn.dealii.org/trunk@17429 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/dofs/dof_renumbering.cc
deal.II/doc/news/changes.h
deal.II/lac/include/lac/sparsity_tools.h
deal.II/lac/source/sparsity_tools.cc

index d1f7ea21d15a6f5a0f212767951d63b7fdd7cf8d..095d619a914852869dc78f9bd11abf07abc1262c 100644 (file)
 //---------------------------------------------------------------------------
 
 
-//TODO:[WB] Unify lots of code of the two Cuthill-McKee dof renumbering functions
-//    This should be rather
-//    straightforward, since all the unified code needs to get is a
-//    sparsity structure, possibly compressed and return a vector
-//    of numbers. Simple task.
-
 #include <base/thread_management.h>
 #include <lac/sparsity_pattern.h>
+#include <lac/sparsity_tools.h>
 #include <lac/compressed_sparsity_pattern.h>
 #include <dofs/dof_accessor.h>
 #include <grid/tria_iterator.h>
@@ -511,179 +506,9 @@ namespace DoFRenumbering
                                     // number was chosen yet
     Assert(new_indices.size() == n_dofs,
           ExcDimensionMismatch(new_indices.size(), n_dofs));
-  
-                                    // store the indices of the dofs renumbered
-                                    // in the last round. Default to starting
-                                    // points
-    std::vector<unsigned int> last_round_dofs (starting_indices);
-  
-                                    // delete disallowed elements
-    for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-      if ((last_round_dofs[i]==DH::invalid_dof_index) ||
-         (last_round_dofs[i]>=n_dofs))
-       last_round_dofs[i] = DH::invalid_dof_index;
-  
-    std::remove_if (last_round_dofs.begin(), last_round_dofs.end(),
-                   std::bind2nd(std::equal_to<unsigned int>(),
-                                DH::invalid_dof_index));
-  
-                                    // now if no valid points remain:
-                                    // find dof with lowest coordination
-                                    // number
-  
-    if (last_round_dofs.size() == 0)
-      {
-       unsigned int starting_point   = DH::invalid_dof_index;
-       unsigned int min_coordination = n_dofs;
-       for (unsigned int row=0; row<n_dofs; ++row) 
-         {
-           unsigned int j;
-
-                                            // loop until we hit the end
-                                            // of this row's entries
-           for (j=sparsity.get_rowstart_indices()[row];
-                j<sparsity.get_rowstart_indices()[row+1]; ++j)
-             if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
-               break;
-                                            // post-condition after loop:
-                                            // coordination, i.e. the number
-                                            // of entries in this row is now
-                                            // j-rowstart[row]
-           if (j-sparsity.get_rowstart_indices()[row] <  min_coordination)
-             {
-               min_coordination = j-sparsity.get_rowstart_indices()[row];
-               starting_point   = row;
-             };
-         };
-      
-                                        // now we still have to care for the
-                                        // case that no dof has a coordination
-                                        // number less than n_dofs. this rather
-                                        // exotic case only happens if we only
-                                        // have one cell, as far as I can see,
-                                        // but there may be others as well.
-                                        //
-                                        // if that should be the case, we can
-                                        // chose an arbitrary dof as starting
-                                        // point, e.g. the one with number zero
-       if (starting_point == DH::invalid_dof_index)
-         starting_point = 0;
-      
-                                        // initialize the first dof
-       last_round_dofs.push_back (starting_point);
-      };
-
 
-                                    // store next free dof index
-    unsigned int next_free_number = 0;
-
-                                    // enumerate the first round dofs
-    for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
-      new_indices[last_round_dofs[i]] = next_free_number++;
-
-    bool all_dofs_renumbered = false;
-
-                                    // now do as many steps as needed to
-                                    // renumber all dofs
-    while (!all_dofs_renumbered) 
-      {
-                                        // store the indices of the dofs to be
-                                        // renumbered in the next round
-       std::vector<unsigned int> next_round_dofs;
-
-                                        // find all neighbors of the
-                                        // dofs numbered in the last
-                                        // round
-       for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-         for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
-              j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
-             break;
-           else
-             next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-      
-                                        // sort dof numbers
-       std::sort (next_round_dofs.begin(), next_round_dofs.end());
-
-                                        // delete multiple entries
-       std::vector<unsigned int>::iterator end_sorted;
-       end_sorted = std::unique (next_round_dofs.begin(), next_round_dofs.end());
-       next_round_dofs.erase (end_sorted, next_round_dofs.end());
-
-                                        // eliminate dofs which are
-                                        // already numbered
-       for (int s=next_round_dofs.size()-1; s>=0; --s)
-         if (new_indices[next_round_dofs[s]] != DH::invalid_dof_index)
-           next_round_dofs.erase (next_round_dofs.begin() + s);
-
-                                        // check whether there are any new
-                                        // dofs
-       all_dofs_renumbered = (next_round_dofs.size() == 0);
-       if (all_dofs_renumbered)
-                                          // end loop if possible
-         continue;
-
-
-                                        // store for each coordination
-                                        // number the dofs with these
-                                        // coordination number
-       std::multimap<unsigned int, int> dofs_by_coordination;
-      
-                                        // find coordination number for
-                                        // each of these dofs
-       for (std::vector<unsigned int>::iterator s=next_round_dofs.begin();
-            s!=next_round_dofs.end(); ++s) 
-         {
-           unsigned int coordination = 0;
-           for (unsigned int j=sparsity.get_rowstart_indices()[*s];
-                j<sparsity.get_rowstart_indices()[*s+1]; ++j)
-             if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
-               break;
-             else
-               ++coordination;
-
-                                            // insert this dof at its
-                                            // coordination number
-           const std::pair<const unsigned int, int> new_entry (coordination, *s);
-           dofs_by_coordination.insert (new_entry);
-         };
-      
-                                        // assign new DoF numbers to
-                                        // the elements of the present
-                                        // front:
-       std::multimap<unsigned int, int>::iterator i;
-       for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
-         new_indices[i->second] = next_free_number++;
-
-                                        // after that: copy this round's
-                                        // dofs for the next round
-       last_round_dofs = next_round_dofs;
-      };
-
-#ifdef DEBUG
-                                    // test for all indices
-                                    // numbered. this mostly tests
-                                    // whether the
-                                    // front-marching-algorithm (which
-                                    // Cuthill-McKee actually is) has
-                                    // reached all points. it should
-                                    // usually do so, but might not for
-                                    // two reasons:
-                                    //
-                                    // - The algorithm above has a bug, or
-                                    // - The domain is not connected and
-                                    // consists of separate parts.
-                                    //
-                                    // In any case, if not all DoFs
-                                    // have been reached, renumbering
-                                    // will not be possible
-    if (std::find (new_indices.begin(), new_indices.end(), DH::invalid_dof_index)
-       !=
-       new_indices.end())
-      Assert (false, ExcRenumberingIncomplete());
-    Assert (next_free_number == n_dofs,
-           ExcRenumberingIncomplete());
-#endif
+    SparsityTools::reorder_Cuthill_McKee (sparsity, new_indices,
+                                         starting_indices);
 
     if (reversed_numbering)
       for (std::vector<unsigned int>::iterator i=new_indices.begin();
@@ -710,145 +535,8 @@ namespace DoFRenumbering
                                     // that no new number was chosen yet
     std::vector<unsigned int> new_indices(n_dofs, DoFHandler<dim>::invalid_dof_index);
   
-                                    // store the indices of the dofs renumbered
-                                    // in the last round. Default to starting
-                                    // points
-    std::vector<unsigned int> last_round_dofs (starting_indices);
-  
-                                    // delete disallowed elements
-    for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-      if ((last_round_dofs[i]==DoFHandler<dim>::invalid_dof_index) ||
-         (last_round_dofs[i]>=n_dofs))
-       last_round_dofs[i] = DoFHandler<dim>::invalid_dof_index;
-  
-    std::remove_if (last_round_dofs.begin(), last_round_dofs.end(),
-                   std::bind2nd(std::equal_to<unsigned int>(),
-                                DoFHandler<dim>::invalid_dof_index));
-  
-                                    // now if no valid points remain:
-                                    // find dof with lowest coordination
-                                    // number
-  
-    if (last_round_dofs.size() == 0)
-      {
-       unsigned int starting_point   = DoFHandler<dim>::invalid_dof_index;
-       unsigned int min_coordination = n_dofs;
-       for (unsigned int row=0; row<n_dofs; ++row) 
-         {
-           unsigned int j;
-           for (j=sparsity.get_rowstart_indices()[row];
-                j<sparsity.get_rowstart_indices()[row+1]; ++j)
-             if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
-               break;
-                                            // post-condition after loop:
-                                            // coordination is now
-                                            // j-rowstart[row]
-           if (j-sparsity.get_rowstart_indices()[row] <  min_coordination)
-             {
-               min_coordination = j-sparsity.get_rowstart_indices()[row];
-               starting_point   = row;
-             };
-         };
-                                        // initialize the first dof
-       last_round_dofs.push_back (starting_point);
-      };
-
-
-                                    // store next free dof index
-    unsigned int next_free_number = 0;
-
-                                    // enumerate the first round dofs
-    for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
-      new_indices[last_round_dofs[i]] = next_free_number++;
-
-    bool all_dofs_renumbered = false;
-
-                                    // now do as many steps as needed to
-                                    // renumber all dofs
-    while (!all_dofs_renumbered) 
-      {
-                                        // store the indices of the dofs to be
-                                        // renumbered in the next round
-       std::vector<unsigned int> next_round_dofs;
-
-                                        // find all neighbors of the
-                                        // dofs numbered in the last
-                                        // round
-       for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-         for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
-              j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
-             break;
-           else
-             next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-      
-                                        // sort dof numbers
-       std::sort (next_round_dofs.begin(), next_round_dofs.end());
-
-                                        // delete multiple entries
-       std::vector<unsigned int>::iterator end_sorted;
-       end_sorted = std::unique (next_round_dofs.begin(), next_round_dofs.end());
-       next_round_dofs.erase (end_sorted, next_round_dofs.end());
-
-                                        // eliminate dofs which are
-                                        // already numbered
-       for (int s=next_round_dofs.size()-1; s>=0; --s)
-         if (new_indices[next_round_dofs[s]] != DoFHandler<dim>::invalid_dof_index)
-           next_round_dofs.erase (next_round_dofs.begin() + s);
-
-                                        // check whether there are any new
-                                        // dofs
-       all_dofs_renumbered = (next_round_dofs.size() == 0);
-       if (all_dofs_renumbered)
-                                          // end loop if possible
-         continue;
-
-
-                                        // store for each coordination
-                                        // number the dofs with these
-                                        // coordination number
-       std::multimap<unsigned int, int> dofs_by_coordination;
-      
-                                        // find coordination number for
-                                        // each of these dofs
-       for (std::vector<unsigned int>::iterator s=next_round_dofs.begin();
-            s!=next_round_dofs.end(); ++s) 
-         {
-           unsigned int coordination = 0;
-           for (unsigned int j=sparsity.get_rowstart_indices()[*s];
-                j<sparsity.get_rowstart_indices()[*s+1]; ++j)
-             if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
-               break;
-             else
-               ++coordination;
-
-                                            // insert this dof at its
-                                            // coordination number
-           const std::pair<const unsigned int, int> new_entry (coordination, *s);
-           dofs_by_coordination.insert (new_entry);
-         };
-      
-                                        ////
-       std::multimap<unsigned int, int>::iterator i;
-       for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
-         new_indices[i->second] = next_free_number++;
-
-                                        // after that: copy this round's
-                                        // dofs for the next round
-       last_round_dofs = next_round_dofs;
-      };
-
-#ifdef DEBUG
-                                    //  test for all indices numbered
-//TODO: Test fails. Do not check before unification.
-    if (std::find (new_indices.begin(), new_indices.end(),
-                  DoFHandler<dim>::invalid_dof_index)
-       !=
-       new_indices.end())
-      Assert (false, ExcRenumberingIncomplete());
-    Assert (next_free_number == n_dofs,
-           ExcRenumberingIncomplete());
-#endif
+    SparsityTools::reorder_Cuthill_McKee (sparsity, new_indices,
+                                         starting_indices);
 
     if (reversed_numbering)
       for (std::vector<unsigned int>::iterator i=new_indices.begin(); i!=new_indices.end(); ++i)
index 1cb171dfedc0b0bd71b01b178b6c5b4417aff54b..d4d4c45c0747bc4cd9d5dcdf6f4be1f5a31637cc 100644 (file)
@@ -301,7 +301,21 @@ inconvenience this causes.
 <ol>
   <li>
   <p>
-  New: The function GridTools::get_face_connectivity_of_cells produces 
+  New: The function SparsityTools::reorder_Cuthill_McKee reorders
+  the nodes of a graph based on their connectivity to other nodes.
+  <br>
+  (WB 2008/10/31)
+  </p>
+
+
+  <li>
+  <p>
+  New: The function GridTools::get_face_connectivity_of_cells produces a
+  sparsity pattern that describes the connectivity of cells of a
+  triangulation based on whether they share common faces.
+  <br>
+  (WB 2008/10/31)
+  </p>
 
   <li>
   <p>
index b448182cefb04f1d86b7437b57990fb1ff0ead61..4fc760555d8a2a96046d18af1fe0925ed82e48b8 100644 (file)
@@ -95,6 +95,81 @@ namespace SparsityTools
                  const unsigned int         n_partitions,
                  std::vector<unsigned int> &partition_indices);
 
+                                  /**
+                                   * For a given sparsity pattern, compute a
+                                   * re-enumeration of row/column indices
+                                   * based on the algorithm by Cuthill-McKee.
+                                   *
+                                   * This algorithm is a graph renumbering
+                                   * algorithm in which we attempt to find a
+                                   * new numbering of all nodes of a graph
+                                   * based on their connectivity to other
+                                   * nodes (i.e. the edges that connect
+                                   * nodes). This connectivity is here
+                                   * represented by the sparsity pattern. In
+                                   * many cases within the library, the nodes
+                                   * represent degrees of freedom and edges
+                                   * are nonzero entries in a matrix,
+                                   * i.e. pairs of degrees of freedom that
+                                   * couple through the action of a bilinear
+                                   * form.
+                                   *
+                                   * The algorithms starts at a node,
+                                   * searches the other nodes for
+                                   * those which are coupled with the one we
+                                   * started with and numbers these in a
+                                   * certain way. It then finds the second
+                                   * level of nodes, namely those that couple
+                                   * with those of the previous level (which
+                                   * were those that coupled with the initial
+                                   * node) and numbers these. And so on. For
+                                   * the details of the algorithm, especially
+                                   * the numbering within each level, we
+                                   * refer the reader to the book of Schwarz
+                                   * (H. R. Schwarz: Methode der finiten
+                                   * Elemente).
+                                   *
+                                   * These algorithms have one major
+                                   * drawback: they require a good starting
+                                   * node, i.e. node that will have number
+                                   * zero in the output array. A starting
+                                   * node forming the initial level of nodes
+                                   * can thus be given by the user, e.g. by
+                                   * exploiting knowledge of the actual
+                                   * topology of the domain. It is also
+                                   * possible to give several starting
+                                   * indices, which may be used to simulate a
+                                   * simple upstream numbering (by giving the
+                                   * inflow nodes as starting values) or to
+                                   * make preconditioning faster (by letting
+                                   * the Dirichlet boundary indices be
+                                   * starting points).
+                                   *
+                                   * If no starting index is given, one is
+                                   * chosen automatically, namely one with
+                                   * the smallest coordination number (the
+                                   * coordination number is the number of
+                                   * other nodes this node couples
+                                   * with). This node is usually located on
+                                   * the boundary of the domain. There is,
+                                   * however, large ambiguity in this when
+                                   * using the hierarchical meshes used in
+                                   * this library, since in most cases the
+                                   * computational domain is not approximated
+                                   * by tilting and deforming elements and by
+                                   * plugging together variable numbers of
+                                   * elements at vertices, but rather by
+                                   * hierarchical refinement. There is
+                                   * therefore a large number of nodes with
+                                   * equal coordination numbers. The
+                                   * renumbering algorithms will therefore
+                                   * not give optimal results.
+                                   */
+  void
+  reorder_Cuthill_McKee (const SparsityPattern     &sparsity,
+                        std::vector<unsigned int> &new_indices,
+                        const std::vector<unsigned int> &starting_indices = std::vector<unsigned int>());
+  
                                   /**
                                    * Exception
                                    */
index 9d844aa8824cf5e569fdcb6be3b39e9c91c381cd..7e8c9b34094492a2e8af1b8cdddf39ff0f0445b0 100644 (file)
@@ -16,6 +16,8 @@
 #include <lac/sparsity_pattern.h>
 #include <lac/sparsity_tools.h>
 
+#include <algorithm>
+
 #ifdef DEAL_II_USE_METIS
 extern "C" {
 #  include <metis.h>
@@ -107,6 +109,195 @@ namespace SparsityTools
 #endif
   }
 
+
+
+  void
+  reorder_Cuthill_McKee (const SparsityPattern     &sparsity,
+                        std::vector<unsigned int> &new_indices,
+                        const std::vector<unsigned int> &starting_indices)
+  {
+    Assert (sparsity.n_rows() == sparsity.n_cols(),
+           ExcDimensionMismatch (sparsity.n_rows(), sparsity.n_cols()));
+    Assert (sparsity.n_rows() == new_indices.size(),
+           ExcDimensionMismatch (sparsity.n_rows(), new_indices.size()));
+    Assert (starting_indices.size() <= sparsity.n_rows(),
+           ExcMessage ("You can't specify more starting indices than there are rows"));
+    for (unsigned int i=0; i<starting_indices.size(); ++i)
+      Assert (starting_indices[i] < sparsity.n_rows(),
+             ExcMessage ("Invalid starting index"));
+    
+                                    // store the indices of the dofs renumbered
+                                    // in the last round. Default to starting
+                                    // points
+    std::vector<unsigned int> last_round_dofs (starting_indices);
+  
+                                    // delete disallowed elements
+    for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+      if ((last_round_dofs[i]==numbers::invalid_unsigned_int) ||
+         (last_round_dofs[i]>=sparsity.n_rows()))
+       last_round_dofs[i] = numbers::invalid_unsigned_int;
+  
+    std::remove_if (last_round_dofs.begin(), last_round_dofs.end(),
+                   std::bind2nd(std::equal_to<unsigned int>(),
+                                numbers::invalid_unsigned_int));
+  
+                                    // now if no valid points remain:
+                                    // find dof with lowest coordination
+                                    // number
+    if (last_round_dofs.size() == 0)
+      {
+       unsigned int starting_point   = numbers::invalid_unsigned_int;
+       unsigned int min_coordination = sparsity.n_rows();
+       for (unsigned int row=0; row<sparsity.n_rows(); ++row) 
+         {
+           unsigned int j;
+
+                                            // loop until we hit the end
+                                            // of this row's entries
+           for (j=sparsity.get_rowstart_indices()[row];
+                j<sparsity.get_rowstart_indices()[row+1]; ++j)
+             if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+               break;
+                                            // post-condition after loop:
+                                            // coordination, i.e. the number
+                                            // of entries in this row is now
+                                            // j-rowstart[row]
+           if (j-sparsity.get_rowstart_indices()[row] <  min_coordination)
+             {
+               min_coordination = j-sparsity.get_rowstart_indices()[row];
+               starting_point   = row;
+             }
+         }
+      
+                                        // now we still have to care for the
+                                        // case that no dof has a coordination
+                                        // number less than sparsity.n_rows(). this rather
+                                        // exotic case only happens if we only
+                                        // have one cell, as far as I can see,
+                                        // but there may be others as well.
+                                        //
+                                        // if that should be the case, we can
+                                        // chose an arbitrary dof as starting
+                                        // point, e.g. the one with number zero
+       if (starting_point == numbers::invalid_unsigned_int)
+         starting_point = 0;
+      
+                                        // initialize the first dof
+       last_round_dofs.push_back (starting_point);
+      }
+
+
+                                    // store next free dof index
+    unsigned int next_free_number = 0;
+
+                                    // enumerate the first round dofs
+    for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
+      new_indices[last_round_dofs[i]] = next_free_number++;
+
+    bool all_dofs_renumbered = false;
+
+                                    // now do as many steps as needed to
+                                    // renumber all dofs
+    while (!all_dofs_renumbered) 
+      {
+                                        // store the indices of the dofs to be
+                                        // renumbered in the next round
+       std::vector<unsigned int> next_round_dofs;
+
+                                        // find all neighbors of the
+                                        // dofs numbered in the last
+                                        // round
+       for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+         for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
+              j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
+           if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+             break;
+           else
+             next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
+      
+                                        // sort dof numbers
+       std::sort (next_round_dofs.begin(), next_round_dofs.end());
+
+                                        // delete multiple entries
+       std::vector<unsigned int>::iterator end_sorted;
+       end_sorted = std::unique (next_round_dofs.begin(), next_round_dofs.end());
+       next_round_dofs.erase (end_sorted, next_round_dofs.end());
+
+                                        // eliminate dofs which are
+                                        // already numbered
+       for (int s=next_round_dofs.size()-1; s>=0; --s)
+         if (new_indices[next_round_dofs[s]] != numbers::invalid_unsigned_int)
+           next_round_dofs.erase (next_round_dofs.begin() + s);
+
+                                        // check whether there are any new
+                                        // dofs
+       all_dofs_renumbered = (next_round_dofs.size() == 0);
+       if (all_dofs_renumbered)
+                                          // end loop if possible
+         continue;
+
+
+                                        // store for each coordination
+                                        // number the dofs with these
+                                        // coordination number
+       std::multimap<unsigned int, int> dofs_by_coordination;
+      
+                                        // find coordination number for
+                                        // each of these dofs
+       for (std::vector<unsigned int>::iterator s=next_round_dofs.begin();
+            s!=next_round_dofs.end(); ++s) 
+         {
+           unsigned int coordination = 0;
+           for (unsigned int j=sparsity.get_rowstart_indices()[*s];
+                j<sparsity.get_rowstart_indices()[*s+1]; ++j)
+             if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+               break;
+             else
+               ++coordination;
+
+                                            // insert this dof at its
+                                            // coordination number
+           const std::pair<const unsigned int, int> new_entry (coordination, *s);
+           dofs_by_coordination.insert (new_entry);
+         }
+      
+                                        // assign new DoF numbers to
+                                        // the elements of the present
+                                        // front:
+       std::multimap<unsigned int, int>::iterator i;
+       for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
+         new_indices[i->second] = next_free_number++;
+
+                                        // after that: copy this round's
+                                        // dofs for the next round
+       last_round_dofs = next_round_dofs;
+      }
+
+                                    // test for all indices
+                                    // numbered. this mostly tests
+                                    // whether the
+                                    // front-marching-algorithm (which
+                                    // Cuthill-McKee actually is) has
+                                    // reached all points. it should
+                                    // usually do so, but might not for
+                                    // two reasons:
+                                    //
+                                    // - The algorithm above has a bug, or
+                                    // - The domain is not connected and
+                                    // consists of separate parts.
+                                    //
+                                    // In any case, if not all DoFs
+                                    // have been reached, renumbering
+                                    // will not be possible
+    Assert ((std::find (new_indices.begin(), new_indices.end(), numbers::invalid_unsigned_int)
+            ==
+            new_indices.end())
+           &&
+           (next_free_number == sparsity.n_rows()),
+           ExcMessage("Some rows have not been renumbered. Maybe "
+                      "the connectivity graph has two or more disconnected "
+                      "subgraphs?"));
+  }
 }
 
 DEAL_II_NAMESPACE_CLOSE

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.