class WrapMGDoFIterator : public T
{
public:
+ WrapMGDoFIterator (const T& t) : T(t) {}
void get_dof_indices (std::vector<unsigned int>& v) const
{
T::get_mg_dof_indices(v);
start = dof_handler.begin_active();
const typename DoFHandler<dim>::cell_iterator
end = dof_handler.end();
-
+
+ unsigned int result =
compute_component_wise<dim,
typename DoFHandler<dim>::active_cell_iterator,
typename DoFHandler<dim>::cell_iterator>(renumbering,
start, end,
component_order_arg);
- if (renumbering.size()!=0)
- dof_handler.renumber_dofs (renumbering);
+ if (result == 0) return;
+
+ Assert (result == dof_handler.n_dofs(),
+ ExcRenumberingIncomplete());
+
+ dof_handler.renumber_dofs (renumbering);
}
template <int dim, class ITERATOR, class ENDITERATOR>
-void
+unsigned int
DoFRenumbering::compute_component_wise (
std::vector<unsigned int>& new_indices,
ITERATOR& start,
if (fe.n_components() == 1)
{
new_indices.resize(0);
- return;
+ return 0;
}
// Copy last argument into a
new_indices[*dof_index] = next_free_index++;
};
-// Assert (next_free_index == dof_handler.n_dofs(),
-// ExcInternalError());
+ return next_free_index;
}
unsigned int level,
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();
+ std::vector<unsigned int> renumbering (dof_handler.n_dofs(level),
+ DoFHandler<dim>::invalid_dof_index);
- // do nothing if the FE has only
- // one component
- if (fe.n_components() == 1)
- return;
+ typedef
+ WrapMGDoFIterator<typename MGDoFHandler<dim>::cell_iterator> ITERATOR;
- std::vector<unsigned int> component_order (component_order_arg);
- if (component_order.size() == 0)
- 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() == fe.n_components(),
- ExcInvalidComponentOrder());
- 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. 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)
- 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 (fe.n_components());
- for (typename MGDoFHandler<dim>::cell_iterator
- cell=dof_handler.begin(level);
- cell!=dof_handler.end(level); ++cell)
- {
- // on each cell: get dof indices
- // and insert them into the global
- // list using their component
- cell->get_mg_dof_indices (local_dof_indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- component_to_dof_map[component_list[i]].push_back (local_dof_indices[i]);
- };
+ typename MGDoFHandler<dim>::cell_iterator
+ istart =dof_handler.begin(level);
+ ITERATOR start = istart;
+ typename MGDoFHandler<dim>::cell_iterator
+ iend = dof_handler.end(level);
+ const ITERATOR end = iend;
+
+ unsigned int result =
+ compute_component_wise<dim, ITERATOR, ITERATOR>(
+ renumbering, start, end, component_order_arg);
+
+ Assert (result == dof_handler.n_dofs(level),
+ ExcRenumberingIncomplete());
- // 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
- //
- // note that we no
- // longer have to care about
- // non-primitive shape functions
- // since the buckets corresponding
- // to the second and following
- // vector components of a
- // non-primitive FE will simply be
- // empty, everything being shoved
- // into the first one
- 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());
- };
-
- // now concatenate all the
- // components in the order the user
- // desired to see
- 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 c=0; c<fe.n_components(); ++c)
- {
- const unsigned int component = component_order[c];
-
- 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;
- dof_index != end_of_component; ++dof_index)
- new_indices[*dof_index] = next_free_index++;
- };
-
- Assert (next_free_index == dof_handler.n_dofs(level),
- ExcInternalError());
-
- dof_handler.renumber_dofs (level, new_indices);
+ if (renumbering.size()!=0)
+ dof_handler.renumber_dofs (renumbering);
}