From: bangerth Date: Fri, 31 Oct 2008 17:24:43 +0000 (+0000) Subject: Split out the function that does the actual Cuthill-McKee algorithm from DoFRenumberi... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b17913a0659b20778fdab65dbffc72520f1f494d;p=dealii-svn.git Split out the function that does the actual Cuthill-McKee algorithm from DoFRenumbering into SparsityTools. git-svn-id: https://svn.dealii.org/trunk@17429 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/dofs/dof_renumbering.cc b/deal.II/deal.II/source/dofs/dof_renumbering.cc index d1f7ea21d1..095d619a91 100644 --- a/deal.II/deal.II/source/dofs/dof_renumbering.cc +++ b/deal.II/deal.II/source/dofs/dof_renumbering.cc @@ -12,14 +12,9 @@ //--------------------------------------------------------------------------- -//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 #include +#include #include #include #include @@ -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 last_round_dofs (starting_indices); - - // delete disallowed elements - for (unsigned int i=0; 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(), - 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 next_round_dofs; - - // find all neighbors of the - // dofs numbered in the last - // round - for (unsigned int i=0; i::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 dofs_by_coordination; - - // find coordination number for - // each of these dofs - for (std::vector::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 new_entry (coordination, *s); - dofs_by_coordination.insert (new_entry); - }; - - // assign new DoF numbers to - // the elements of the present - // front: - std::multimap::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::iterator i=new_indices.begin(); @@ -710,145 +535,8 @@ namespace DoFRenumbering // that no new number was chosen yet std::vector new_indices(n_dofs, DoFHandler::invalid_dof_index); - // store the indices of the dofs renumbered - // in the last round. Default to starting - // points - std::vector last_round_dofs (starting_indices); - - // delete disallowed elements - for (unsigned int i=0; i::invalid_dof_index) || - (last_round_dofs[i]>=n_dofs)) - last_round_dofs[i] = DoFHandler::invalid_dof_index; - - std::remove_if (last_round_dofs.begin(), last_round_dofs.end(), - std::bind2nd(std::equal_to(), - DoFHandler::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::invalid_dof_index; - unsigned int min_coordination = n_dofs; - for (unsigned int row=0; row next_round_dofs; - - // find all neighbors of the - // dofs numbered in the last - // round - for (unsigned int i=0; i::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::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 dofs_by_coordination; - - // find coordination number for - // each of these dofs - for (std::vector::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 new_entry (coordination, *s); - dofs_by_coordination.insert (new_entry); - }; - - //// - std::multimap::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::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::iterator i=new_indices.begin(); i!=new_indices.end(); ++i) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 1cb171dfed..d4d4c45c07 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -301,7 +301,21 @@ inconvenience this causes.
  1. - 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. +
    + (WB 2008/10/31) +

    + + +
  2. +

    + 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. +
    + (WB 2008/10/31) +

  3. diff --git a/deal.II/lac/include/lac/sparsity_tools.h b/deal.II/lac/include/lac/sparsity_tools.h index b448182cef..4fc760555d 100644 --- a/deal.II/lac/include/lac/sparsity_tools.h +++ b/deal.II/lac/include/lac/sparsity_tools.h @@ -95,6 +95,81 @@ namespace SparsityTools const unsigned int n_partitions, std::vector &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 &new_indices, + const std::vector &starting_indices = std::vector()); + /** * Exception */ diff --git a/deal.II/lac/source/sparsity_tools.cc b/deal.II/lac/source/sparsity_tools.cc index 9d844aa882..7e8c9b3409 100644 --- a/deal.II/lac/source/sparsity_tools.cc +++ b/deal.II/lac/source/sparsity_tools.cc @@ -16,6 +16,8 @@ #include #include +#include + #ifdef DEAL_II_USE_METIS extern "C" { # include @@ -107,6 +109,195 @@ namespace SparsityTools #endif } + + + void + reorder_Cuthill_McKee (const SparsityPattern &sparsity, + std::vector &new_indices, + const std::vector &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 last_round_dofs (starting_indices); + + // delete disallowed elements + for (unsigned int i=0; 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(), + 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 next_round_dofs; + + // find all neighbors of the + // dofs numbered in the last + // round + for (unsigned int i=0; i::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 dofs_by_coordination; + + // find coordination number for + // each of these dofs + for (std::vector::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 new_entry (coordination, *s); + dofs_by_coordination.insert (new_entry); + } + + // assign new DoF numbers to + // the elements of the present + // front: + std::multimap::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