std::vector<SparsityPattern::size_type> &new_indices,
const std::vector<SparsityPattern::size_type> &starting_indices = std::vector<SparsityPattern::size_type>()) DEAL_II_DEPRECATED;
+ /**
+ * For a given sparsity pattern, compute a re-enumeration of row/column
+ * indices in a hierarchical way, similar to what
+ * DoFRenumbering::hierarchical does for degrees of freedom on
+ * hierarchically refined meshes.
+ *
+ * This algorithm first selects a node with the minimum number of neighbors
+ * and puts that node and its direct neighbors into one chunk. Next, it
+ * selects one of the neighbors of the already selected nodes, adds the node
+ * and its direct neighbors that are not part of one of the previous chunks,
+ * into the next. After this sweep, neighboring nodes are grouped
+ * together. To ensure a similar grouping on a more global level, this
+ * grouping is called recursively on the groups so formed. The recursion
+ * stops when no further grouping is possible. Eventually, the ordering
+ * obtained by this method passes through the indices represented in the
+ * sparsity pattern in a z-like way.
+ *
+ * If the graph has two or more unconnected components, the algorithm will
+ * number each component consecutively, starting with the components with
+ * the lowest number of nodes.
+ */
+ void
+ reorder_hierarchical (const DynamicSparsityPattern &sparsity,
+ std::vector<DynamicSparsityPattern::size_type> &new_indices);
#ifdef DEAL_II_WITH_MPI
/**
#include <algorithm>
#include <functional>
+#include <set>
#ifdef DEAL_II_WITH_MPI
#include <deal.II/base/utilities.h>
for (unsigned int row=0; row<sparsity.n_rows(); ++row)
{
for (SparsityPattern::iterator it=sparsity.begin(row); it!=sparsity.end(row)
- && it->is_valid_entry() ; ++it)
+ && it->is_valid_entry() ; ++it)
dsp.add(row, it->column());
}
reorder_Cuthill_McKee(dsp, new_indices, starting_indices);
+ namespace internal
+ {
+ void
+ reorder_hierarchical (const DynamicSparsityPattern &connectivity,
+ std::vector<DynamicSparsityPattern::size_type> &renumbering)
+ {
+ AssertDimension (connectivity.n_rows(), connectivity.n_cols());
+ AssertDimension (connectivity.n_rows(), renumbering.size());
+
+ std::vector<types::global_dof_index> touched_nodes(connectivity.n_rows(),
+ numbers::invalid_dof_index);
+ std::set<types::global_dof_index> current_neighbors;
+ std::vector<std::vector<types::global_dof_index> > groups;
+
+ // This outer loop is typically traversed only once, unless the global
+ // graph is not connected
+ while (true)
+ {
+ // Find cell with the minimal number of neighbors (typically a corner
+ // node when based on FEM meshes). If no cell is left, we are done.
+ std::pair<types::global_dof_index,types::global_dof_index> min_neighbors
+ (numbers::invalid_unsigned_int, numbers::invalid_unsigned_int);
+ for (types::global_dof_index i=0; i<touched_nodes.size(); ++i)
+ if (touched_nodes[i] == numbers::invalid_unsigned_int)
+ if (connectivity.row_length(i) < min_neighbors.second)
+ min_neighbors = std::make_pair(i, connectivity.row_length(i));
+ if (min_neighbors.first == numbers::invalid_unsigned_int)
+ break;
+
+ current_neighbors.clear();
+ current_neighbors.insert(min_neighbors.first);
+ while (!current_neighbors.empty())
+ {
+ // Find cell with minimum number of untouched neighbors among the
+ // next set of possible neighbors
+ min_neighbors = std::make_pair (numbers::invalid_unsigned_int,
+ numbers::invalid_unsigned_int);
+ for (std::set<types::global_dof_index>::iterator it=current_neighbors.begin();
+ it != current_neighbors.end(); ++it)
+ {
+ AssertThrow(touched_nodes[*it] == numbers::invalid_unsigned_int,
+ ExcInternalError());
+ types::global_dof_index active_row_length = 0;
+ for (CompressedSimpleSparsityPattern::row_iterator rowit
+ = connectivity.row_begin(*it);
+ rowit != connectivity.row_end(*it); ++rowit)
+ if (touched_nodes[*rowit] == numbers::invalid_unsigned_int)
+ ++active_row_length;
+ if (active_row_length < min_neighbors.second)
+ min_neighbors = std::make_pair(*it, active_row_length);
+ }
+ // Among the set of cells with the minimal number of neighbors,
+ // choose the one with the largest number of touched neighbors,
+ // i.e., the one with the largest row length
+ const types::global_dof_index best_row_length = min_neighbors.second;
+ for (std::set<types::global_dof_index>::iterator it=current_neighbors.begin();
+ it != current_neighbors.end(); ++it)
+ {
+ types::global_dof_index active_row_length = 0;
+ for (CompressedSimpleSparsityPattern::row_iterator rowit
+ = connectivity.row_begin(*it);
+ rowit != connectivity.row_end(*it); ++rowit)
+ if (touched_nodes[*rowit] == numbers::invalid_unsigned_int)
+ ++active_row_length;
+ if (active_row_length == best_row_length)
+ if (connectivity.row_length(*it) > min_neighbors.second)
+ min_neighbors = std::make_pair(*it, connectivity.row_length(*it));
+ }
+
+ // Add the pivot cell to the current list
+ groups.push_back(std::vector<types::global_dof_index>());
+ std::vector<types::global_dof_index> &new_entries = groups.back();
+ new_entries.push_back(min_neighbors.first);
+ touched_nodes[min_neighbors.first] = groups.size()-1;
+
+ // Add all direct neighbors of the pivot cell not yet touched to
+ // the current list
+ for (CompressedSimpleSparsityPattern::row_iterator it
+ = connectivity.row_begin(min_neighbors.first);
+ it != connectivity.row_end(min_neighbors.first); ++it)
+ {
+ if (touched_nodes[*it] == numbers::invalid_unsigned_int)
+ {
+ new_entries.push_back(*it);
+ touched_nodes[*it] = groups.size()-1;
+ }
+ }
+
+ // Add all neighbors of the current list not yet touched to the
+ // set of possible next pivots. Delete the entries of the current
+ // list from the set of possible next pivots.
+ for (types::global_dof_index i=0; i<new_entries.size(); ++i)
+ {
+ for (CompressedSimpleSparsityPattern::row_iterator it
+ = connectivity.row_begin(new_entries[i]);
+ it != connectivity.row_end(new_entries[i]); ++it)
+ if (touched_nodes[*it] == numbers::invalid_unsigned_int)
+ current_neighbors.insert(*it);
+ current_neighbors.erase(new_entries[i]);
+ }
+ }
+ }
+ // If the number of groups is smaller than the number of nodes, we can
+ // continue by recursively calling this method
+ if (groups.size() < connectivity.n_rows())
+ {
+ // Form the connectivity of the groups
+ DynamicSparsityPattern connectivity_next(groups.size(),
+ groups.size());
+ for (types::global_dof_index i=0; i<groups.size(); ++i)
+ for (types::global_dof_index col=0; col<groups[i].size(); ++col)
+ for (CompressedSimpleSparsityPattern::row_iterator it
+ = connectivity.row_begin(groups[i][col]);
+ it != connectivity.row_end(groups[i][col]); ++it)
+ {
+ if (touched_nodes[*it] != i)
+ connectivity_next.add(i, touched_nodes[*it]);
+ }
+
+ // Recursively call the reordering
+ std::vector<types::global_dof_index> renumbering_next(groups.size());
+ reorder_hierarchical(connectivity_next, renumbering_next);
+
+ // Renumber the indices group by group according to the incoming
+ // ordering for the groups
+ for (types::global_dof_index i=0,count=0; i<groups.size(); ++i)
+ for (types::global_dof_index col=0; col<groups[renumbering_next[i]].size(); ++col, ++count)
+ renumbering[count] = groups[renumbering_next[i]][col];
+ }
+ else
+ {
+ // All groups should have size one and no more recursion is possible,
+ // so use the numbering of the groups
+ for (types::global_dof_index i=0,count=0; i<groups.size(); ++i)
+ for (types::global_dof_index col=0; col<groups[i].size(); ++col, ++count)
+ renumbering[count] = groups[i][col];
+ }
+ }
+ }
+
+ void
+ reorder_hierarchical (const DynamicSparsityPattern &connectivity,
+ std::vector<DynamicSparsityPattern::size_type> &renumbering)
+ {
+ // the internal renumbering keeps the numbering the wrong way around (but
+ // we cannot invert the numbering inside that method because it is used
+ // recursively), so invert it here
+ internal::reorder_hierarchical(connectivity, renumbering);
+ renumbering = Utilities::invert_permutation(renumbering);
+ }
+
+
+
#ifdef DEAL_II_WITH_MPI
template <class DSP_t>
void distribute_sparsity_pattern(DSP_t &dsp,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2008 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// apply SparsityTools::reorder_hierarchical to a graph that consists
+// of two non-connected parts.
+
+#include "../tests.h"
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+
+#include <iomanip>
+#include <fstream>
+
+int main ()
+{
+ std::ofstream logfile("output");
+ logfile.setf(std::ios::fixed);
+ deallog << std::setprecision(3);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ DynamicSparsityPattern dsp (8,8);
+ for (unsigned int i=0; i<8; ++i)
+ dsp.add (i,i);
+
+ // create a graph with components 0,2 and 1,3 that are disconnected
+ dsp.add (0,2);
+ dsp.add (2,0);
+
+ dsp.add (1,3);
+ dsp.add (3,1);
+
+ // couple all indices between 3 and 7 except 4
+ for (unsigned int i=3; i<7; ++i)
+ for (unsigned int j=3; j<7; ++j)
+ if (i != 4 && j != 4)
+ dsp.add(i,j);
+
+ // indices 4,7 couple
+ dsp.add (4,7);
+ dsp.add (7,4);
+
+ // now find permutation
+ std::vector<types::global_dof_index> permutation(8);
+ SparsityTools::reorder_hierarchical (dsp, permutation);
+
+ for (unsigned int i=0; i<permutation.size(); ++i)
+ deallog << permutation[i] << std::endl;
+}
+
+
+