* may be difficult, however, and in many cases will not justify the effort.
*
*
- * <h3>Component-wise numbering</h3>
+ * <h3>Component-wise and block-wise numberings</h3>
*
* 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.
+ * different components. For example, one may want to sort degree of freedom for
+ * a Stokes discretization so that we first get all velocities and then all
+ * the pressures so that the resulting matrix naturally decomposes into a
+ * $2\times 2$ system.
*
* 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
+ * the order of indices within each component, 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
*
* The component_wise() function allows not only to honor enumeration based on
* vector components, but also allows to group together vector components into
- * "blocks". See @ref GlossComponent vs @ref GlossBlock for a description of
- * the difference.
+ * "blocks" using a defaulted argument to the various DoFRenumber::component_wise()
+ * functinos (see @ref GlossComponent vs @ref GlossBlock for a description of
+ * the difference). The blocks designated through this argument may, but do not
+ * have to be, equal to the blocks that the finite element reports. For example,
+ * a typical Stokes element would be
+ * @code
+ * FESystem<dim> stokes_fe (FE_Q<dim>(2), dim, // dim velocities
+ * FE_Q<dim>(1), 1); // one pressure
+ * @endcode
+ * This element has <code>dim+1</code> vector components and equally many
+ * blocks. However, one may want to consider the velocities as one logical
+ * block so that all velocity degrees of freedom are enumerated the same
+ * way, independent of whether they are $x$- or $y$-velocities. This is done,
+ * for example, in step-20 and step-22 as well as several other tutorial programs.
+ *
+ * On the other hand, if you really want to use block structure reported
+ * by the finite element itself (a case that is often the case if you have
+ * finite elements that have multiple vector components, e.g. the FE_RaviartThomas
+ * or FE_Nedelec elements) then you can use the DoFRenumber::block_wise instead
+ * of the DoFRenumbering::component_wise functions.
*
*
* <h3>Cell-wise numbering</h3>
* target component for dofs with
* component @p i in the
* FESystem. Naming the same
- * component more than once is
+ * target component more than once is
* possible and results in a
* blocking of several components
* into one. This is discussed in
template <int dim>
void
component_wise (MGDoFHandler<dim>& dof_handler,
- unsigned int level,
+ const unsigned int level,
const std::vector<unsigned int>& target_component = std::vector<unsigned int>());
* @}
*/
+ /**
+ * @name Block-wise numberings
+ * @{
+ */
+
+ /**
+ * Sort the degrees of freedom by
+ * vector block. The
+ * numbering within each
+ * block is not touched, so a
+ * degree of freedom with index
+ * $i$, belonging to some
+ * block, and another degree
+ * of freedom with index $j$
+ * belonging to the same
+ * block 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$.
+ */
+ template <int dim, int spacedim>
+ void
+ block_wise (DoFHandler<dim,spacedim> &dof_handler);
+
+
+ /**
+ * Sort the degrees of freedom by
+ * block. It does the same
+ * thing as the above function.
+ */
+ template <int dim>
+ void
+ block_wise (hp::DoFHandler<dim> &dof_handler);
+
+ /**
+ * Sort the degrees of freedom by
+ * block. It does the same
+ * thing as the above function,
+ * only that it does this for one
+ * single level of a multi-level
+ * discretization. The
+ * non-multigrid part of the
+ * MGDoFHandler is not touched.
+ */
+ template <int dim>
+ void
+ block_wise (MGDoFHandler<dim> &dof_handler,
+ const unsigned int level);
+
+
+ /**
+ * Sort the degrees of freedom by
+ * block. It does the same
+ * thing as the previous
+ * functions, but more: it
+ * renumbers not only every level
+ * of the multigrid part, but
+ * also the global,
+ * i.e. non-multigrid components.
+ */
+ template <int dim>
+ void
+ block_wise (MGDoFHandler<dim> &dof_handler);
+
+ /**
+ * Computes the renumbering
+ * vector needed by the
+ * block_wise()
+ * functions. Does not perform
+ * the renumbering on the
+ * DoFHandler dofs but returns
+ * the renumbering vector.
+ */
+ template <int dim, int spacedim, class ITERATOR, class ENDITERATOR>
+ unsigned int
+ compute_block_wise (std::vector<unsigned int>& new_dof_indices,
+ const ITERATOR& start,
+ const ENDITERATOR& end);
+
+ /**
+ * @}
+ */
+
/**
* @name Various cell-wise numberings
* @{
}
// now we've got all indices sorted
- // into buckets labelled with their
+ // into buckets labeled by their
// target component number. we've
// only got to traverse this list
// and assign the new indices
return next_free_index;
}
+
+
+ template <int dim, int spacedim>
+ void
+ block_wise (DoFHandler<dim,spacedim> &dof_handler)
+ {
+ std::vector<unsigned int> renumbering (dof_handler.n_locally_owned_dofs(),
+ DoFHandler<dim>::invalid_dof_index);
+
+ typedef
+ internal::WrapDoFIterator<typename DoFHandler<dim,spacedim>
+ ::active_cell_iterator>
+ ITERATOR;
+
+ typename DoFHandler<dim,spacedim>::active_cell_iterator
+ istart = dof_handler.begin_active();
+ ITERATOR start = istart;
+ const typename DoFHandler<dim,spacedim>::cell_iterator
+ end = dof_handler.end();
+
+ const unsigned int result =
+ compute_block_wise<dim, spacedim, ITERATOR,
+ typename DoFHandler<dim,spacedim>::cell_iterator>
+ (renumbering, start, end);
+ if (result == 0)
+ return;
+
+ // verify that the last numbered
+ // degree of freedom is either
+ // equal to the number of degrees
+ // of freedom in total (the
+ // sequential case) or in the
+ // distributed case at least
+ // makes sense
+ Assert ((result == dof_handler.n_locally_owned_dofs())
+ ||
+ ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs())
+ &&
+ (result <= dof_handler.n_dofs())),
+ ExcRenumberingIncomplete());
+
+ dof_handler.renumber_dofs (renumbering);
+ }
+
+
+
+ template <int dim>
+ void
+ block_wise (hp::DoFHandler<dim> &dof_handler)
+ {
+ std::vector<unsigned int> renumbering (dof_handler.n_dofs(),
+ hp::DoFHandler<dim>::invalid_dof_index);
+
+ typedef
+ internal::WrapDoFIterator<typename hp::DoFHandler<dim>::active_cell_iterator> ITERATOR;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ istart = dof_handler.begin_active();
+ ITERATOR start = istart;
+ const typename hp::DoFHandler<dim>::cell_iterator
+ end = dof_handler.end();
+
+ const unsigned int result =
+ compute_block_wise<dim, dim, ITERATOR,
+ typename hp::DoFHandler<dim>::cell_iterator>(renumbering,
+ start, end);
+
+ if (result == 0)
+ return;
+
+ Assert (result == dof_handler.n_dofs(),
+ ExcRenumberingIncomplete());
+
+ dof_handler.renumber_dofs (renumbering);
+ }
+
+
+
+ template <int dim>
+ void
+ block_wise (MGDoFHandler<dim> &dof_handler,
+ const unsigned int level)
+ {
+ std::vector<unsigned int> renumbering (dof_handler.n_dofs(level),
+ DoFHandler<dim>::invalid_dof_index);
+
+ typedef
+ internal::WrapMGDoFIterator<typename MGDoFHandler<dim>::cell_iterator> ITERATOR;
+
+ 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;
+
+ const unsigned int result =
+ compute_block_wise<dim, dim, ITERATOR, ITERATOR>(
+ renumbering, start, end);
+
+ if (result == 0) return;
+
+ Assert (result == dof_handler.n_dofs(level),
+ ExcRenumberingIncomplete());
+
+ if (renumbering.size()!=0)
+ dof_handler.renumber_dofs (level, renumbering);
+ }
+
+
+
+ template <int dim>
+ void
+ block_wise (MGDoFHandler<dim> &dof_handler)
+ {
+ // renumber the non-MG part of
+ // the DoFHandler in parallel to
+ // the MG part. Because
+ // MGDoFHandler::renumber_dofs
+ // uses the user flags we can't
+ // run renumbering on individual
+ // levels in parallel to the
+ // other levels
+ void (*non_mg_part) (DoFHandler<dim> &)
+ = &block_wise<dim>;
+ Threads::Task<>
+ task = Threads::new_task (non_mg_part, dof_handler);
+
+ for (unsigned int level=0; level<dof_handler.get_tria().n_levels(); ++level)
+ block_wise (dof_handler, level);
+
+ task.join();
+ }
+
+
+
+ template <int dim, int spacedim, class ITERATOR, class ENDITERATOR>
+ unsigned int
+ compute_block_wise (std::vector<unsigned int>& new_indices,
+ const ITERATOR & start,
+ const ENDITERATOR& end)
+ {
+ const hp::FECollection<dim,spacedim>
+ fe_collection (start->get_dof_handler().get_fe ());
+
+ // do nothing if the FE has only
+ // one component
+ if (fe_collection.n_blocks() == 1)
+ {
+ new_indices.resize(0);
+ return 0;
+ }
+
+ // vector to hold the dof indices on
+ // the cell we visit at a time
+ std::vector<unsigned int> local_dof_indices;
+
+ // prebuilt list to which block
+ // a given dof on a cell
+ // should go.
+ std::vector<std::vector<unsigned int> > block_list (fe_collection.size());
+ for (unsigned int f=0; f<fe_collection.size(); ++f)
+ {
+ const FiniteElement<dim,spacedim> & fe = fe_collection[f];
+ block_list[f].resize(fe.dofs_per_cell);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ block_list[f][i]
+ = fe.system_to_block_index(i).first;
+ }
+
+ // set up a map where for each
+ // block the respective degrees
+ // of freedom are collected.
+ //
+ // note that this map is sorted by
+ // block but that within each
+ // block 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> >
+ block_to_dof_map (fe_collection.n_blocks());
+ for (ITERATOR cell=start; cell!=end; ++cell)
+ if (cell->is_locally_owned())
+ {
+ // on each cell: get dof indices
+ // and insert them into the global
+ // list using their component
+ const unsigned int fe_index = cell->active_fe_index();
+ const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell;
+ local_dof_indices.resize (dofs_per_cell);
+ cell.get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (start->get_dof_handler().locally_owned_dofs().is_element(local_dof_indices[i]))
+ block_to_dof_map[block_list[fe_index][i]].
+ push_back (local_dof_indices[i]);
+ }
+
+ // now we've got all indices sorted
+ // into buckets labeled by their
+ // target block 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 block=0; block<fe_collection.n_blocks();
+ ++block)
+ {
+ std::sort (block_to_dof_map[block].begin(),
+ block_to_dof_map[block].end());
+ block_to_dof_map[block]
+ .erase (std::unique (block_to_dof_map[block].begin(),
+ block_to_dof_map[block].end()),
+ block_to_dof_map[block].end());
+ }
+
+ // calculate the number of locally owned
+ // DoFs per bucket
+ const unsigned int n_buckets = fe_collection.n_blocks();
+ std::vector<unsigned int> shifts(n_buckets);
+
+ if (const parallel::distributed::Triangulation<dim,spacedim> * tria
+ = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
+ (&start->get_dof_handler().get_tria())))
+ {
+#ifdef DEAL_II_USE_P4EST
+ std::vector<unsigned int> local_dof_count(n_buckets);
+
+ for (unsigned int c=0; c<n_buckets; ++c)
+ local_dof_count[c] = block_to_dof_map[c].size();
+
+
+ // gather information from all CPUs
+ std::vector<unsigned int>
+ all_dof_counts(fe_collection.n_components() *
+ Utilities::MPI::n_mpi_processes (tria->get_communicator()));
+
+ MPI_Allgather ( &local_dof_count[0], n_buckets, MPI_UNSIGNED, &all_dof_counts[0],
+ n_buckets, MPI_UNSIGNED, tria->get_communicator());
+
+ for (unsigned int i=0; i<n_buckets; ++i)
+ Assert (all_dof_counts[n_buckets*tria->locally_owned_subdomain()+i]
+ ==
+ local_dof_count[i],
+ ExcInternalError());
+
+ //calculate shifts
+ unsigned int cumulated = 0;
+ for (unsigned int c=0; c<n_buckets; ++c)
+ {
+ shifts[c]=cumulated;
+ for (types::subdomain_id i=0; i<tria->locally_owned_subdomain(); ++i)
+ shifts[c] += all_dof_counts[c+n_buckets*i];
+ for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes (tria->get_communicator()); ++i)
+ cumulated += all_dof_counts[c+n_buckets*i];
+ }
+#else
+ (void)tria;
+ Assert (false, ExcInternalError());
+#endif
+ }
+ else
+ {
+ shifts[0] = 0;
+ for (unsigned int c=1; c<fe_collection.n_blocks(); ++c)
+ shifts[c] = shifts[c-1] + block_to_dof_map[c-1].size();
+ }
+
+
+
+
+ // now concatenate all the
+ // components in the order the user
+ // desired to see
+ unsigned int next_free_index = 0;
+ for (unsigned int block=0; block<fe_collection.n_blocks(); ++block)
+ {
+ const typename std::vector<unsigned int>::const_iterator
+ begin_of_component = block_to_dof_map[block].begin(),
+ end_of_component = block_to_dof_map[block].end();
+
+ next_free_index = shifts[block];
+
+ for (typename std::vector<unsigned int>::const_iterator
+ dof_index = begin_of_component;
+ dof_index != end_of_component; ++dof_index)
+ {
+ Assert (start->get_dof_handler().locally_owned_dofs()
+ .index_within_set(*dof_index)
+ <
+ new_indices.size(),
+ ExcInternalError());
+ new_indices[start->get_dof_handler().locally_owned_dofs()
+ .index_within_set(*dof_index)]
+ = next_free_index++;
+ }
+ }
+
+ return next_free_index;
+ }
+
+
+
namespace
{
// helper function for hierarchical()
projects / dealii-svn.git / commitdiff