//---------------------------------------------------------------------------
-//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>
// 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();
// 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)
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
*/
#include <lac/sparsity_pattern.h>
#include <lac/sparsity_tools.h>
+#include <algorithm>
+
#ifdef DEAL_II_USE_METIS
extern "C" {
# include <metis.h>
#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