/**
* Sort the degrees of freedom by
- * component. The numbering within
+ * vector component. The numbering within
* each component is not touched,
* so a degree of freedom with index
* $i$, belonging to some component,
* in the finite element, and has to
* contain all numbers counted from
* zero onwards. If
- * you ommit this argument, the same
+ * you omit 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.
+ * If one of the base finite
+ * elements from which the global
+ * finite element under
+ * consideration here, is a
+ * non-primitive one, i.e. its
+ * shape functions have more than
+ * one non-zero component, then
+ * it is not possible to
+ * associate these degrees of
+ * freedom with a single vector
+ * component. In this case, they
+ * are associated with the first
+ * vector component to which they
+ * belong.
+ *
+ * For finite elements with only
+ * one component, or a single
+ * non-primitive base element,
+ * this function is the identity
+ * operation.
*/
template <int dim>
static void
/**
* 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 @p{(0, 3, 2, 1)}, then all
- * indices of component @p{0} will be
- * before those of component @p{3}, before
- * those of component @p{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.
+ * component. It does the same
+ * thing as the above function,
+ * only that it does this for one
+ * level of a multi-level
+ * discretization.
*/
template <int dim>
static void
template <int dim>
-void DoFRenumbering::component_wise (DoFHandler<dim> &dof_handler,
- const std::vector<unsigned int> &component_order_arg)
+void
+DoFRenumbering::component_wise (DoFHandler<dim> &dof_handler,
+ const std::vector<unsigned int> &component_order_arg)
{
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ const FiniteElement<dim> &fe = dof_handler.get_fe();
// do nothing if the FE has only
// one component
- if (dof_handler.get_fe().n_components() == 1)
+ if (fe.n_components() == 1)
return;
std::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)
+ for (unsigned int i=0; i<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(),
+ Assert (component_order.size() == fe.n_components(),
ExcInvalidComponentOrder());
- for (unsigned int i=0; i<dof_handler.get_fe().n_components(); ++i)
+ for (unsigned int i=0; i<fe.n_components(); ++i)
Assert (std::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
std::vector<unsigned int> local_dof_indices(dofs_per_cell);
- // prebuilt list to which component
- // a given dof on a cell belongs
+
+ // prebuilt list to which component
+ // a given dof on a cell
+ // belongs. note that we get into
+ // trouble here if the shape
+ // function is not primitive, since
+ // then there is no single vector
+ // component to which it
+ // belongs. in this case, assign it
+ // to the first vector component to
+ // which it belongs
std::vector<unsigned int> component_list (dofs_per_cell);
for (unsigned int i=0; i<component_list.size(); ++i)
- component_list[i] = dof_handler.get_fe().system_to_component_index(i).first;
+ if (fe.is_primitive(i))
+ component_list[i] = fe.system_to_component_index(i).first;
+ else
+ {
+ const unsigned int base_element =
+ fe.system_to_base_index(i).first.first;
+ const unsigned int base_multiplicity =
+ fe.system_to_base_index(i).first.second;
+ // sum up the number of
+ // components all the
+ // elements before that have
+ unsigned int c=0;
+ for (unsigned int b=0; b<base_element; ++b)
+ c += fe.base_element(b).n_components() *
+ fe.element_multiplicity(b);
+ for (unsigned int m=0; m<base_multiplicity; ++m)
+ c += fe.base_element(base_element).n_components();
+ // then associate this degree
+ // of freedom with this
+ // component
+ component_list[i] = c;
+ };
- // set up a map where for each component
- // the respective degrees of freedom are
- // collected.
+ // 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
- std::vector<std::vector<unsigned int> > component_to_dof_map (dof_handler.get_fe().n_components());
- for (typename DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active();
+ // 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
+ std::vector<std::vector<unsigned int> >
+ component_to_dof_map (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
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
+ // 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
+ // 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)
+ for (unsigned int component=0; component<fe.n_components(); ++component)
{
std::sort (component_to_dof_map[component].begin(),
component_to_dof_map[component].end());
- component_to_dof_map[component].erase (std::unique (component_to_dof_map[component].begin(),
- component_to_dof_map[component].end()),
- component_to_dof_map[component].end());
+ component_to_dof_map[component]
+ .erase (std::unique (component_to_dof_map[component].begin(),
+ component_to_dof_map[component].end()),
+ component_to_dof_map[component].end());
};
unsigned int next_free_index = 0;
std::vector<unsigned int> new_indices (dof_handler.n_dofs(),
DoFHandler<dim>::invalid_dof_index);
- for (unsigned int component=0; component<dof_handler.get_fe().n_components(); ++component)
+ for (unsigned int component=0; component<fe.n_components(); ++component)
{
const typename std::vector<unsigned int>::const_iterator
begin_of_component = component_to_dof_map[component].begin(),
end_of_component = component_to_dof_map[component].end();
- for (typename std::vector<unsigned int>::const_iterator dof_index = begin_of_component;
+ for (typename std::vector<unsigned int>::const_iterator
+ dof_index = begin_of_component;
dof_index != end_of_component; ++dof_index)
new_indices[*dof_index] = next_free_index++;
};
const std::vector<unsigned int> &component_order_arg)
{
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ const FiniteElement<dim> &fe = dof_handler.get_fe();
// do nothing if the FE has only
// one component
- if (dof_handler.get_fe().n_components() == 1)
+ if (fe.n_components() == 1)
return;
std::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)
+ for (unsigned int i=0; i<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(),
+ Assert (component_order.size() == fe.n_components(),
ExcInvalidComponentOrder());
- for (unsigned int i=0; i<dof_handler.get_fe().n_components(); ++i)
+ for (unsigned int i=0; i<fe.n_components(); ++i)
Assert (std::find (component_order.begin(), component_order.end(), i)
!= component_order.end(),
ExcInvalidComponentOrder ());
std::vector<unsigned int> local_dof_indices(dofs_per_cell);
// prebuilt list to which component
// a given dof on a cell belongs
+ // prebuilt list to which component
+ // a given dof on a cell
+ // belongs. note that we get into
+ // trouble here if the shape
+ // function is not primitive, since
+ // then there is no single vector
+ // component to which it
+ // belongs. in this case, assign it
+ // to the first vector component to
+ // which it belongs
std::vector<unsigned int> component_list (dofs_per_cell);
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.
+ if (fe.is_primitive(i))
+ component_list[i] = fe.system_to_component_index(i).first;
+ else
+ {
+ const unsigned int base_element =
+ fe.system_to_base_index(i).first.first;
+ const unsigned int base_multiplicity =
+ fe.system_to_base_index(i).first.second;
+ // sum up the number of
+ // components all the
+ // elements before that have
+ unsigned int c=0;
+ for (unsigned int b=0; b<base_element; ++b)
+ c += fe.base_element(b).n_components() *
+ fe.element_multiplicity(b);
+ for (unsigned int m=0; m<base_multiplicity; ++m)
+ c += fe.base_element(base_element).n_components();
+ // then associate this degree
+ // of freedom with this
+ // component
+ component_list[i] = c;
+ };
+
+ // 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
- std::vector<std::vector<unsigned int> > component_to_dof_map (dof_handler.get_fe().n_components());
+ // 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
+ std::vector<std::vector<unsigned int> > component_to_dof_map (fe.n_components());
for (typename MGDoFHandler<dim>::cell_iterator cell=dof_handler.begin(level);
cell!=dof_handler.end(level); ++cell)
{
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
+ // 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
+ // 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)
+ for (unsigned int component=0; component<fe.n_components(); ++component)
{
std::sort (component_to_dof_map[component].begin(),
component_to_dof_map[component].end());
unsigned int next_free_index = 0;
std::vector<unsigned int> new_indices (dof_handler.n_dofs(level),
DoFHandler<dim>::invalid_dof_index);
- for (unsigned int component=0; component<dof_handler.get_fe().n_components(); ++component)
+ for (unsigned int component=0; component<fe.n_components(); ++component)
{
const typename std::vector<unsigned int>::const_iterator
begin_of_component = component_to_dof_map[component].begin(),