]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add DoFRenumbering::block_wise functions.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 20 Sep 2012 13:37:48 +0000 (13:37 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 20 Sep 2012 13:37:48 +0000 (13:37 +0000)
git-svn-id: https://svn.dealii.org/trunk@26565 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/dofs/dof_renumbering.h
deal.II/source/dofs/dof_renumbering.cc
deal.II/source/dofs/dof_renumbering.inst.in

index bb85d16439b07ab6e4a1af5722e16ddbca16144c..6b7d838b1f07e598671b71105722b42efc20fb6b 100644 (file)
@@ -144,18 +144,21 @@ DEAL_II_NAMESPACE_OPEN
  * 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
@@ -163,8 +166,26 @@ DEAL_II_NAMESPACE_OPEN
  *
  * 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>
@@ -700,7 +721,7 @@ namespace DoFRenumbering
                                     * 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
@@ -760,7 +781,7 @@ namespace DoFRenumbering
   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>());
 
 
@@ -799,6 +820,89 @@ namespace DoFRenumbering
    * @}
    */
 
+  /**
+   * @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
    * @{
index 99bec3b82325ce10f5ca1b0845ff9161c0cc9b0e..01eaf7e9f312fc5c4e67603579f5f242d2d7b4ab 100644 (file)
@@ -781,7 +781,7 @@ namespace DoFRenumbering
         }
 
                                      // 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
@@ -899,6 +899,315 @@ namespace DoFRenumbering
     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()
index 769af6b944c1748ada8b79b8ebbfc1f44257b78c..fd9feb3cd4603a910a5907d62076c61c87f2e606 100644 (file)
@@ -44,6 +44,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
          (DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
           const std::vector<unsigned int>&);
 
+       template
+         void block_wise<deal_II_dimension,deal_II_space_dimension>
+         (DoFHandler<deal_II_dimension,deal_II_space_dimension>&);
+
        template
          void subdomain_wise<DoFHandler<deal_II_dimension,deal_II_space_dimension> >
          (DoFHandler<deal_II_dimension,deal_II_space_dimension> &);
@@ -169,6 +173,19 @@ namespace DoFRenumbering
       (MGDoFHandler<deal_II_dimension>&,
        const std::vector<unsigned int>&);
 
+    template
+      void block_wise<deal_II_dimension>
+      (hp::DoFHandler<deal_II_dimension>&);
+
+    template
+      void block_wise<deal_II_dimension>
+      (MGDoFHandler<deal_II_dimension>&,
+       unsigned int);
+
+    template
+      void block_wise<deal_II_dimension>
+      (MGDoFHandler<deal_II_dimension>&);
+
     template
       void hierarchical<deal_II_dimension>
       (DoFHandler<deal_II_dimension>&);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.