void make_sparsity_pattern (const unsigned int level,
dSMatrixStruct &sparsity) const;
+ /**
+ * Make up the constraint matrix which
+ * is used to condensate the
+ * system matrices and to prolong
+ * the solution vectors from the true
+ * degrees of freedom also to the
+ * constraint nodes. The constraints
+ * apply only to the given level.
+ *
+ * Since this method does not make sense in
+ * one dimension, the functions returns
+ * immediately after clearing the
+ * constraint matrix.
+ * For more than one dimension, the matrix
+ * is cleared before usage. The constraint
+ * matrix is closed anyway, no matter of the
+ * dimension.
+ *
+ * To condense a given sparsity pattern,
+ * use #ConstraintMatrix::condense#.
+ *
+ * This function uses the user flags for
+ * the faces. If you need the user flags,
+ * store them beforehand.
+ */
+ void make_constraint_matrix (const unsigned int level,
+ ConstraintMatrix &cm) const;
+
/*--------------------------------------*/
/**
* set all indices to #-1#.
*/
void init (const unsigned int coarsest_level,
- const unsigned int n_levels,
+ const unsigned int finest_level,
const unsigned int dofs_per_vertex);
/**
int get_index (const unsigned int level,
const unsigned int dof_number,
const unsigned int dofs_per_vertex) const;
+
+ /**
+ * Return the index of the coarsest
+ * level this vertex lives on.
+ */
+ unsigned int get_coarsest_level () const;
+
+ /**
+ * Return the index of the finest
+ * level this vertex lives on.
+ */
+ unsigned int get_finest_level () const;
/**
* Exception.
*/
DeclException1 (ExcInvalidLevel,
int,
- << "The given level number " << arg1 << " is below the "
- << "coarsest level this vertex lives on.");
+ << "The given level number " << arg1 << " is outside "
+ << "the range of levels this vertex lives on.");
/**
* Exception.
*/
*/
unsigned int coarsest_level;
+ /**
+ * Finest level this level lives on.
+ * This is mostly used for error
+ * checking but can also be accessed
+ * by the function #get_finest_level#.
+ */
+ unsigned int finest_level;
+
/**
* Array holding the indices.
*/
int *indices;
};
-
- /**
- * Reserve enough space for the MG dof
- * indices for a given triangulation.
- */
- void reserve_space ();
/**
* Distribute dofs on the given cell,
*/
unsigned int distribute_dofs_on_cell (active_cell_iterator &cell,
unsigned int next_free_dof);
+
+ /**
+ * Actually do the renumbering prepared
+ * by the #renumber_dofs# function on
+ * the given #level#. Since
+ * this is dimension specific, we
+ * need to have another function.
+ *
+ * #new_numbers# is an array of integers
+ * with size equal to the number of dofs
+ * on the present level. It stores the new
+ * indices after renumbering in the
+ * order of the old indices.
+ */
+ void do_renumbering (const unsigned int level,
+ const vector<int> &new_numbers);
+
+ /**
+ * Reserve enough space for the MG dof
+ * indices for a given triangulation.
+ */
+ void reserve_space ();
/**
* Space to store the DoF numbers for the
const unsigned int dof_number,
const unsigned int dofs_per_vertex,
const unsigned int index) {
- Assert (level >= coarsest_level, ExcInvalidLevel(level));
+ Assert ((level >= coarsest_level) && (level <= finest_level),
+ ExcInvalidLevel(level));
Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ());
indices[(level-coarsest_level)*dofs_per_vertex + dof_number] = index;
int MGDoFHandler<dim>::MGVertexDoFs::get_index (const unsigned int level,
const unsigned int dof_number,
const unsigned int dofs_per_vertex) const {
- Assert (level >= coarsest_level, ExcInvalidLevel(level));
+ Assert ((level >= coarsest_level) && (level <= finest_level),
+ ExcInvalidLevel(level));
Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ());
return indices[(level-coarsest_level)*dofs_per_vertex + dof_number];
void make_sparsity_pattern (const unsigned int level,
dSMatrixStruct &sparsity) const;
+ /**
+ * Make up the constraint matrix which
+ * is used to condensate the
+ * system matrices and to prolong
+ * the solution vectors from the true
+ * degrees of freedom also to the
+ * constraint nodes. The constraints
+ * apply only to the given level.
+ *
+ * Since this method does not make sense in
+ * one dimension, the functions returns
+ * immediately after clearing the
+ * constraint matrix.
+ * For more than one dimension, the matrix
+ * is cleared before usage. The constraint
+ * matrix is closed anyway, no matter of the
+ * dimension.
+ *
+ * To condense a given sparsity pattern,
+ * use #ConstraintMatrix::condense#.
+ *
+ * This function uses the user flags for
+ * the faces. If you need the user flags,
+ * store them beforehand.
+ */
+ void make_constraint_matrix (const unsigned int level,
+ ConstraintMatrix &cm) const;
+
/*--------------------------------------*/
/**
* set all indices to #-1#.
*/
void init (const unsigned int coarsest_level,
- const unsigned int n_levels,
+ const unsigned int finest_level,
const unsigned int dofs_per_vertex);
/**
int get_index (const unsigned int level,
const unsigned int dof_number,
const unsigned int dofs_per_vertex) const;
+
+ /**
+ * Return the index of the coarsest
+ * level this vertex lives on.
+ */
+ unsigned int get_coarsest_level () const;
+
+ /**
+ * Return the index of the finest
+ * level this vertex lives on.
+ */
+ unsigned int get_finest_level () const;
/**
* Exception.
*/
DeclException1 (ExcInvalidLevel,
int,
- << "The given level number " << arg1 << " is below the "
- << "coarsest level this vertex lives on.");
+ << "The given level number " << arg1 << " is outside "
+ << "the range of levels this vertex lives on.");
/**
* Exception.
*/
*/
unsigned int coarsest_level;
+ /**
+ * Finest level this level lives on.
+ * This is mostly used for error
+ * checking but can also be accessed
+ * by the function #get_finest_level#.
+ */
+ unsigned int finest_level;
+
/**
* Array holding the indices.
*/
int *indices;
};
-
- /**
- * Reserve enough space for the MG dof
- * indices for a given triangulation.
- */
- void reserve_space ();
/**
* Distribute dofs on the given cell,
*/
unsigned int distribute_dofs_on_cell (active_cell_iterator &cell,
unsigned int next_free_dof);
+
+ /**
+ * Actually do the renumbering prepared
+ * by the #renumber_dofs# function on
+ * the given #level#. Since
+ * this is dimension specific, we
+ * need to have another function.
+ *
+ * #new_numbers# is an array of integers
+ * with size equal to the number of dofs
+ * on the present level. It stores the new
+ * indices after renumbering in the
+ * order of the old indices.
+ */
+ void do_renumbering (const unsigned int level,
+ const vector<int> &new_numbers);
+
+ /**
+ * Reserve enough space for the MG dof
+ * indices for a given triangulation.
+ */
+ void reserve_space ();
/**
* Space to store the DoF numbers for the
const unsigned int dof_number,
const unsigned int dofs_per_vertex,
const unsigned int index) {
- Assert (level >= coarsest_level, ExcInvalidLevel(level));
+ Assert ((level >= coarsest_level) && (level <= finest_level),
+ ExcInvalidLevel(level));
Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ());
indices[(level-coarsest_level)*dofs_per_vertex + dof_number] = index;
int MGDoFHandler<dim>::MGVertexDoFs::get_index (const unsigned int level,
const unsigned int dof_number,
const unsigned int dofs_per_vertex) const {
- Assert (level >= coarsest_level, ExcInvalidLevel(level));
+ Assert ((level >= coarsest_level) && (level <= finest_level),
+ ExcInvalidLevel(level));
Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ());
return indices[(level-coarsest_level)*dofs_per_vertex + dof_number];
#include <fe/fe.h>
#include <lac/dsmatrix.h>
+#include <algorithm>
+
+
/* ------------------------ MGVertexDoFs ----------------------------------- */
template <int dim>
MGDoFHandler<dim>::MGVertexDoFs::MGVertexDoFs () :
coarsest_level (1<<30),
+ finest_level (0),
indices (0)
{};
template <int dim>
void MGDoFHandler<dim>::MGVertexDoFs::init (const unsigned int cl,
- const unsigned int n_levels,
+ const unsigned int fl,
const unsigned int dofs_per_vertex) {
Assert (indices == 0, ExcInternalError());
coarsest_level = cl;
+ finest_level = fl;
+
+ const unsigned int n_levels = finest_level-coarsest_level + 1;
indices = new int[n_levels * dofs_per_vertex];
Assert (indices != 0, ExcNoMemory ());
+template <int dim>
+inline
+unsigned int MGDoFHandler<dim>::MGVertexDoFs::get_coarsest_level () const {
+ return coarsest_level;
+};
+
+
+
+template <int dim>
+inline
+unsigned int MGDoFHandler<dim>::MGVertexDoFs::get_finest_level () const {
+ return finest_level;
+};
+
+
+template <int dim>
+void MGDoFHandler<dim>::renumber_dofs (const unsigned int level,
+ const RenumberingMethod method,
+ const bool use_constraints,
+ const vector<int> &starting_points) {
+ // make the connection graph
+ dSMatrixStruct sparsity (n_dofs(level), max_couplings_between_dofs());
+ make_sparsity_pattern (level, sparsity);
+
+ if (use_constraints)
+ {
+ ConstraintMatrix constraints;
+ make_constraint_matrix (level, constraints);
+ constraints.condense (sparsity);
+ };
+
+ int n_dofs = sparsity.n_rows();
+ // store the new dof numbers; -1 means
+ // that no new number was chosen yet
+ //
+ // the commented line is what would be the
+ // correct way to do, but gcc2.8 chokes
+ // over that. The other lines are a
+ // workaround
+// vector<int> new_number(n_dofs, -1);
+ vector<int> new_number;
+ new_number.resize (n_dofs, -1);
+
+ // store the indices of the dofs renumbered
+ // in the last round. Default to starting
+ // points
+ vector<int> last_round_dofs (starting_points);
+
+ // delete disallowed elements
+ for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+ if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
+ last_round_dofs[i] = -1;
+
+ remove_if (last_round_dofs.begin(), last_round_dofs.end(),
+ bind2nd(equal_to<int>(), -1));
+
+ // now if no valid points remain:
+ // find dof with lowest coordination
+ // number
+
+ if (last_round_dofs.size() == 0)
+ {
+ int starting_point = -1;
+ unsigned int min_coordination = n_dofs;
+ for (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] == -1)
+ 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
+ int next_free_number = 0;
+
+ // enumerate the first round dofs
+ for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
+ new_number[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
+ vector<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] == -1)
+ break;
+ else
+ next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
+
+ // sort dof numbers
+ sort (next_round_dofs.begin(), next_round_dofs.end());
+
+ // delete multiple entries
+ vector<int>::iterator end_sorted;
+ end_sorted = 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_number[next_round_dofs[s]] != -1)
+ next_round_dofs.erase (&next_round_dofs[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
+ multimap<unsigned int, int> dofs_by_coordination;
+
+ // find coordination number for
+ // each of these dofs
+ for (vector<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] == -1)
+ break;
+ else
+ ++coordination;
+
+ // insert this dof at its
+ // coordination number
+ const pair<const unsigned int, int> new_entry (coordination, *s);
+ dofs_by_coordination.insert (new_entry);
+ };
+
+ ////
+ multimap<unsigned int, int>::iterator i;
+ for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i)
+ new_number[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
+ if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
+ Assert (false, ExcRenumberingIncomplete());
+ Assert (next_free_number == n_dofs,
+ ExcRenumberingIncomplete());
+#endif
+
+ switch (method)
+ {
+ case Cuthill_McKee:
+ break;
+ case reverse_Cuthill_McKee:
+ {
+ for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
+ *i = n_dofs-*i;
+ break;
+ };
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ // actually perform renumbering;
+ // this is dimension specific and
+ // thus needs an own function
+ do_renumbering (level, new_number);
+};
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+void MGDoFHandler<1>::do_renumbering (const unsigned int level,
+ const vector<int> &new_numbers) {
+ Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
+
+ // note that we can not use cell iterators
+ // in this function since then we would
+ // renumber the dofs on the interface of
+ // two cells more than once. Anyway, this
+ // ways it's not only more correct but also
+ // faster
+ for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
+ i!=mg_vertex_dofs.end(); ++i)
+ // if the present vertex lives on
+ // the present level
+ if ((i->get_coarsest_level() <= level) &&
+ (i->get_finest_level() >= level))
+ for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+ i->set_index (level, d, selected_fe->dofs_per_vertex,
+ new_numbers[i->get_index (level, d,
+ selected_fe->dofs_per_vertex)]);
+
+ for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
+ i!=mg_levels[level]->line_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+void MGDoFHandler<2>::do_renumbering (const unsigned int level,
+ const vector<int> &new_numbers) {
+ Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
+
+ for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
+ i!=mg_vertex_dofs.end(); ++i)
+ // if the present vertex lives on
+ // the present level
+ if ((i->get_coarsest_level() <= level) &&
+ (i->get_finest_level() >= level))
+ for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+ i->set_index (level, d, selected_fe->dofs_per_vertex,
+ new_numbers[i->get_index (level, d,
+ selected_fe->dofs_per_vertex)]);
+
+ for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
+ i!=mg_levels[level]->line_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+ for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
+ i!=mg_levels[level]->quad_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+};
+
+#endif
+
+
+
#if deal_II_dimension == 1
Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
mg_vertex_dofs[vertex].init (min_level[vertex],
- max_level[vertex]-min_level[vertex]+1,
+ max_level[vertex],
selected_fe->dofs_per_vertex);
};
};
Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
mg_vertex_dofs[vertex].init (min_level[vertex],
- max_level[vertex]-min_level[vertex]+1,
+ max_level[vertex],
selected_fe->dofs_per_vertex);
};
};
#include <fe/fe.h>
#include <lac/dsmatrix.h>
+#include <algorithm>
+
+
/* ------------------------ MGVertexDoFs ----------------------------------- */
template <int dim>
MGDoFHandler<dim>::MGVertexDoFs::MGVertexDoFs () :
coarsest_level (1<<30),
+ finest_level (0),
indices (0)
{};
template <int dim>
void MGDoFHandler<dim>::MGVertexDoFs::init (const unsigned int cl,
- const unsigned int n_levels,
+ const unsigned int fl,
const unsigned int dofs_per_vertex) {
Assert (indices == 0, ExcInternalError());
coarsest_level = cl;
+ finest_level = fl;
+
+ const unsigned int n_levels = finest_level-coarsest_level + 1;
indices = new int[n_levels * dofs_per_vertex];
Assert (indices != 0, ExcNoMemory ());
+template <int dim>
+inline
+unsigned int MGDoFHandler<dim>::MGVertexDoFs::get_coarsest_level () const {
+ return coarsest_level;
+};
+
+
+
+template <int dim>
+inline
+unsigned int MGDoFHandler<dim>::MGVertexDoFs::get_finest_level () const {
+ return finest_level;
+};
+
+
+template <int dim>
+void MGDoFHandler<dim>::renumber_dofs (const unsigned int level,
+ const RenumberingMethod method,
+ const bool use_constraints,
+ const vector<int> &starting_points) {
+ // make the connection graph
+ dSMatrixStruct sparsity (n_dofs(level), max_couplings_between_dofs());
+ make_sparsity_pattern (level, sparsity);
+
+ if (use_constraints)
+ {
+ ConstraintMatrix constraints;
+ make_constraint_matrix (level, constraints);
+ constraints.condense (sparsity);
+ };
+
+ int n_dofs = sparsity.n_rows();
+ // store the new dof numbers; -1 means
+ // that no new number was chosen yet
+ //
+ // the commented line is what would be the
+ // correct way to do, but gcc2.8 chokes
+ // over that. The other lines are a
+ // workaround
+// vector<int> new_number(n_dofs, -1);
+ vector<int> new_number;
+ new_number.resize (n_dofs, -1);
+
+ // store the indices of the dofs renumbered
+ // in the last round. Default to starting
+ // points
+ vector<int> last_round_dofs (starting_points);
+
+ // delete disallowed elements
+ for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+ if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
+ last_round_dofs[i] = -1;
+
+ remove_if (last_round_dofs.begin(), last_round_dofs.end(),
+ bind2nd(equal_to<int>(), -1));
+
+ // now if no valid points remain:
+ // find dof with lowest coordination
+ // number
+
+ if (last_round_dofs.size() == 0)
+ {
+ int starting_point = -1;
+ unsigned int min_coordination = n_dofs;
+ for (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] == -1)
+ 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
+ int next_free_number = 0;
+
+ // enumerate the first round dofs
+ for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
+ new_number[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
+ vector<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] == -1)
+ break;
+ else
+ next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
+
+ // sort dof numbers
+ sort (next_round_dofs.begin(), next_round_dofs.end());
+
+ // delete multiple entries
+ vector<int>::iterator end_sorted;
+ end_sorted = 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_number[next_round_dofs[s]] != -1)
+ next_round_dofs.erase (&next_round_dofs[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
+ multimap<unsigned int, int> dofs_by_coordination;
+
+ // find coordination number for
+ // each of these dofs
+ for (vector<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] == -1)
+ break;
+ else
+ ++coordination;
+
+ // insert this dof at its
+ // coordination number
+ const pair<const unsigned int, int> new_entry (coordination, *s);
+ dofs_by_coordination.insert (new_entry);
+ };
+
+ ////
+ multimap<unsigned int, int>::iterator i;
+ for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i)
+ new_number[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
+ if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
+ Assert (false, ExcRenumberingIncomplete());
+ Assert (next_free_number == n_dofs,
+ ExcRenumberingIncomplete());
+#endif
+
+ switch (method)
+ {
+ case Cuthill_McKee:
+ break;
+ case reverse_Cuthill_McKee:
+ {
+ for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
+ *i = n_dofs-*i;
+ break;
+ };
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ // actually perform renumbering;
+ // this is dimension specific and
+ // thus needs an own function
+ do_renumbering (level, new_number);
+};
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+void MGDoFHandler<1>::do_renumbering (const unsigned int level,
+ const vector<int> &new_numbers) {
+ Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
+
+ // note that we can not use cell iterators
+ // in this function since then we would
+ // renumber the dofs on the interface of
+ // two cells more than once. Anyway, this
+ // ways it's not only more correct but also
+ // faster
+ for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
+ i!=mg_vertex_dofs.end(); ++i)
+ // if the present vertex lives on
+ // the present level
+ if ((i->get_coarsest_level() <= level) &&
+ (i->get_finest_level() >= level))
+ for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+ i->set_index (level, d, selected_fe->dofs_per_vertex,
+ new_numbers[i->get_index (level, d,
+ selected_fe->dofs_per_vertex)]);
+
+ for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
+ i!=mg_levels[level]->line_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+void MGDoFHandler<2>::do_renumbering (const unsigned int level,
+ const vector<int> &new_numbers) {
+ Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
+
+ for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
+ i!=mg_vertex_dofs.end(); ++i)
+ // if the present vertex lives on
+ // the present level
+ if ((i->get_coarsest_level() <= level) &&
+ (i->get_finest_level() >= level))
+ for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+ i->set_index (level, d, selected_fe->dofs_per_vertex,
+ new_numbers[i->get_index (level, d,
+ selected_fe->dofs_per_vertex)]);
+
+ for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
+ i!=mg_levels[level]->line_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+ for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
+ i!=mg_levels[level]->quad_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+};
+
+#endif
+
+
+
#if deal_II_dimension == 1
Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
mg_vertex_dofs[vertex].init (min_level[vertex],
- max_level[vertex]-min_level[vertex]+1,
+ max_level[vertex],
selected_fe->dofs_per_vertex);
};
};
Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
mg_vertex_dofs[vertex].init (min_level[vertex],
- max_level[vertex]-min_level[vertex]+1,
+ max_level[vertex],
selected_fe->dofs_per_vertex);
};
};