* may be difficult, however, and in many cases will not justify the effort.
*
*
- * \subsection{Multigrid DoF numbering}
+ * \section{Componentwise numbering}
+ *
+ * For finite elements composed of several base elements using the #FESystem#
+ * class, or for elements which provide several components themselves, it
+ * may be of interest to sort the DoF indices by component. This will then
+ * bring out the block matrix structure, since otherwise the degrees of freedom
+ * are numbered cell-wise without taking into account that they may belong to
+ * different components.
+ *
+ * This kind of numbering may be obtained by calling the #component_wise# function
+ * of this class. Since it does not touch the order of indices within each, it
+ * may be worthwhile to first renumber using the Cuthill-McKee or a similar
+ * algorithm and afterwards renumbering component-wise. This will bring out the
+ * matrix structure and additionally have a good numbering within each block.
+ *
+ *
+ * \section{Multigrid DoF numbering}
*
* Most algorithms also work on multigrid degree of freedom numberings. Refer
* to the actual function declarations to get more information on this.
const bool reversed_numbering = false,
const vector<int> &starting_indices = vector<int> ());
+ /**
+ * Sort the degrees of freedom by
+ * component. The numbering within
+ * each component is not touched,
+ * so a degree of freedom with index
+ * $i$, belonging to some component,
+ * and another degree of freedom
+ * with index $j$ belonging to the same
+ * component will be assigned new
+ * indices $n(i)$ and $n(j)$ with
+ * $n(i)<n(j)$ if $i<j$ and
+ * $n(i)>n(j)$ if $i>j$.
+ *
+ * You may want to give the order in
+ * which the components are to be ordered
+ * (e.g. if the second argument contains
+ * the numbers #(0, 3, 2, 1)#, then all
+ * indices of component #0# will be
+ * before those of component #3#, before
+ * those of component #2#, ...). The
+ * length of this list has to be the
+ * same as the number of components
+ * in the finite element, and has to
+ * contain all numbers counted from
+ * zero onwards. If
+ * you ommit this argument, the same
+ * order as given by the finite element
+ * is used.
+ *
+ * For finite elements with only one
+ * component, this function is the
+ * identity operation.
+ */
+ template <int dim>
+ static void component_wise (DoFHandler<dim> &dof_handler,
+ const vector<unsigned int> &component_order = vector<unsigned int>());
+
/**
* Exception
*/
DeclException0 (ExcRenumberingIncomplete);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidComponentOrder);
};
/* $Id$ */
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
#include <grid/dof.h>
#include <grid/mg_dof.h>
#include <grid/dof_constraints.h>
+#include <fe/fe.h>
#include <numerics/dof_renumbering.h>
#include <lac/sparsematrix.h>
+template <int dim>
+void DoFRenumbering::component_wise (DoFHandler<dim> &dof_handler,
+ const vector<unsigned int> &component_order_arg)
+{
+ const unsigned int total_dofs = dof_handler.get_fe().total_dofs;
+
+ // do nothing if the FE has only
+ // one component
+ if (dof_handler.get_fe().n_components == 1)
+ return;
+
+ vector<unsigned int> component_order (component_order_arg);
+ if (component_order.size() == 0)
+ for (unsigned int i=0; i<dof_handler.get_fe().n_components; ++i)
+ component_order.push_back (i);
+
+ // check whether the component list has
+ // the right length and contains all
+ // component numbers
+ Assert (component_order.size() == dof_handler.get_fe().n_components,
+ ExcInvalidComponentOrder());
+ for (unsigned int i=0; i<dof_handler.get_fe().n_components; ++i)
+ Assert (find (component_order.begin(), component_order.end(), i)
+ != component_order.end(),
+ ExcInvalidComponentOrder ());
+
+ // vector to hold the dof indices on
+ // the cell we visit at a time
+ vector<int> local_dof_indices(total_dofs);
+ // prebuilt list to which component
+ // a given dof on a cell belongs
+ vector<unsigned int> component_list (total_dofs);
+ for (unsigned int i=0; i<component_list.size(); ++i)
+ component_list[i] = dof_handler.get_fe().system_to_component_index(i).first;
+
+ // set up a map where for each component
+ // the respective degrees of freedom are
+ // collected.
+ //
+ // note that this map is sorted by component
+ // but that within each component it is NOT
+ // sorted by dof index. note also that some
+ // dof indices are entered multiply, so we
+ // will have to take care of that
+ vector<vector<int> > component_to_dof_map (dof_handler.get_fe().n_components);
+ for (typename DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active();
+ cell!=dof_handler.end(); ++cell)
+ {
+ // on each cell: get dof indices
+ // and insert them into the global
+ // list using their component
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<total_dofs; ++i)
+ component_to_dof_map[component_list[i]].push_back (local_dof_indices[i]);
+ };
+
+ // now we've got all indices sorted into
+ // buckets labelled with their component
+ // number. we've only got to traverse this
+ // list and assign the new indices
+ //
+ // however, we first want to sort the
+ // indices entered into the buckets to
+ // preserve the order within each component
+ // and during this also remove duplicate
+ // entries
+ for (unsigned int component=0; component<dof_handler.get_fe().n_components; ++component)
+ {
+ sort (component_to_dof_map[component].begin(),
+ component_to_dof_map[component].end());
+ component_to_dof_map[component].erase (unique (component_to_dof_map[component].begin(),
+ component_to_dof_map[component].end()),
+ component_to_dof_map[component].end());
+ };
+
+ int next_free_index = 0;
+ vector<int> new_indices (dof_handler.n_dofs(), -1);
+ for (unsigned int component=0; component<dof_handler.get_fe().n_components; ++component)
+ {
+ const typename vector<int>::const_iterator
+ begin_of_component = component_to_dof_map[component].begin(),
+ end_of_component = component_to_dof_map[component].end();
+
+ for (typename vector<int>::const_iterator dof_index = begin_of_component;
+ dof_index != end_of_component; ++dof_index)
+ new_indices[*dof_index] = next_free_index++;
+ };
+
+ Assert (next_free_index == static_cast<int>(dof_handler.n_dofs()),
+ ExcInternalError());
+
+ dof_handler.renumber_dofs (new_indices);
+};
+
+
+
const bool reversed_numbering,
const vector<int> &starting_indices);
+template
+void DoFRenumbering::component_wise (DoFHandler<deal_II_dimension> &dof_handler,
+ const vector<unsigned int> &component_order_arg);