From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Tue, 22 Sep 2015 16:17:07 +0000 (+0200)
Subject: Enable MGTransferPrebuilt with parallel::distributed::Vector and Trilinos.
X-Git-Tag: v8.4.0-rc2~378^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3d2632b247c342a44f857a3bcc48fad1f7a180a2;p=dealii.git

Enable MGTransferPrebuilt with parallel::distributed::Vector and Trilinos.
Fix MGTransfer for 64 bit indices (and computations beyond 32 bit).
---

diff --git a/cmake/config/template-arguments.in b/cmake/config/template-arguments.in
index ee4444296b..a28f358e92 100644
--- a/cmake/config/template-arguments.in
+++ b/cmake/config/template-arguments.in
@@ -57,6 +57,8 @@ VECTORS_WITH_MATRIX := { Vector<double>;
                     BlockVector<float>;
                     BlockVector<long double>;
 
+                    parallel::distributed::Vector<double>;
+
                     @DEAL_II_EXPAND_TRILINOS_VECTOR@;
                     @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
                   }
diff --git a/include/deal.II/multigrid/mg_transfer.h b/include/deal.II/multigrid/mg_transfer.h
index 03fa244a43..13add34b41 100644
--- a/include/deal.II/multigrid/mg_transfer.h
+++ b/include/deal.II/multigrid/mg_transfer.h
@@ -23,6 +23,7 @@
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/block_sparsity_pattern.h>
 #include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/parallel_vector.h>
 
 #include <deal.II/lac/vector_memory.h>
 
@@ -58,6 +59,21 @@ namespace internal
   };
 
 #ifdef DEAL_II_WITH_TRILINOS
+  template <typename Number>
+  struct MatrixSelector<parallel::distributed::Vector<Number> >
+  {
+    typedef ::dealii::TrilinosWrappers::SparsityPattern Sparsity;
+    typedef ::dealii::TrilinosWrappers::SparseMatrix Matrix;
+
+    template <class DSP, class DH>
+    static void reinit(Matrix &matrix, Sparsity &, int level, const DSP &dsp, DH &dh)
+    {
+      matrix.reinit(dh.locally_owned_mg_dofs(level+1),
+                    dh.locally_owned_mg_dofs(level),
+                    dsp, MPI_COMM_WORLD, true);
+    }
+
+  };
   template <>
   struct MatrixSelector<dealii::TrilinosWrappers::MPI::Vector>
   {
@@ -256,7 +272,7 @@ private:
    * The data is organized as follows: one vector per level. Each element of
    * these vectors contains first the global index, then the level index.
    */
-  std::vector<std::vector<std::pair<types::global_dof_index, unsigned int> > >
+  std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > >
   copy_indices;
 
   /**
@@ -266,7 +282,7 @@ private:
    *
    * Organization of the data is like for @p copy_indices_mine.
    */
-  std::vector<std::vector<std::pair<types::global_dof_index, unsigned int> > >
+  std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > >
   copy_indices_to_me;
 
   /**
@@ -276,7 +292,7 @@ private:
    *
    * Organization of the data is like for @p copy_indices_mine.
    */
-  std::vector<std::vector<std::pair<types::global_dof_index, unsigned int> > >
+  std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > >
   copy_indices_from_me;
 
 
diff --git a/include/deal.II/multigrid/mg_transfer.templates.h b/include/deal.II/multigrid/mg_transfer.templates.h
index b1ff432a3b..1af86a22af 100644
--- a/include/deal.II/multigrid/mg_transfer.templates.h
+++ b/include/deal.II/multigrid/mg_transfer.templates.h
@@ -96,6 +96,35 @@ namespace
       }
   }
 
+  /**
+   * Adjust vectors on all levels to correct size.  Here, we just count the
+   * numbers of degrees of freedom on each level and @p reinit each level
+   * vector to this length.
+   */
+  template <int dim, typename number, int spacedim>
+  void
+  reinit_vector (const dealii::DoFHandler<dim,spacedim> &mg_dof,
+                 const std::vector<unsigned int> &,
+                 MGLevelObject<parallel::distributed::Vector<number> > &v)
+  {
+    const parallel::Triangulation<dim,spacedim> *tria =
+      (dynamic_cast<const parallel::Triangulation<dim,spacedim>*>
+       (&mg_dof.get_tria()));
+
+    for (unsigned int level=v.min_level();
+         level<=v.max_level(); ++level)
+      {
+        const IndexSet vector_index_set = v[level].locally_owned_elements();
+        if (vector_index_set.size() != mg_dof.locally_owned_mg_dofs(level).size() ||
+            mg_dof.locally_owned_mg_dofs(level) != vector_index_set)
+          v[level].reinit(mg_dof.locally_owned_mg_dofs(level), tria != 0 ?
+                          tria->get_communicator() : MPI_COMM_SELF);
+        else
+          v[level] = 0.;
+      }
+  }
+
+
 #ifdef DEAL_II_WITH_TRILINOS
   /**
    * Adjust vectors on all levels to correct size.  Here, we just count the
@@ -148,7 +177,7 @@ MGTransferPrebuilt<VECTOR>::copy_to_mg (
       --level;
       VECTOR &dst_level = dst[level];
 
-      typedef std::vector<std::pair<types::global_dof_index, unsigned int> >::const_iterator IT;
+      typedef std::vector<std::pair<types::global_dof_index, types::global_dof_index> >::const_iterator IT;
       for (IT i= copy_indices[level].begin();
            i != copy_indices[level].end(); ++i)
         dst_level(i->second) = src(i->first);
@@ -186,7 +215,7 @@ MGTransferPrebuilt<VECTOR>::copy_from_mg(
   dst = 0;
   for (unsigned int level=0; level<mg_dof_handler.get_tria().n_global_levels(); ++level)
     {
-      typedef std::vector<std::pair<types::global_dof_index, unsigned int> >::const_iterator IT;
+      typedef std::vector<std::pair<types::global_dof_index, types::global_dof_index> >::const_iterator IT;
 
       // First copy all indices local to this process
       if (constraints==0)
@@ -230,7 +259,7 @@ MGTransferPrebuilt<VECTOR>::copy_from_mg_add (
   // functions
   for (unsigned int level=0; level<mg_dof_handler.get_tria().n_global_levels(); ++level)
     {
-      typedef std::vector<std::pair<types::global_dof_index, unsigned int> >::const_iterator IT;
+      typedef std::vector<std::pair<types::global_dof_index, types::global_dof_index> >::const_iterator IT;
       if (constraints==0)
         for (IT i= copy_indices[level].begin();
              i != copy_indices[level].end(); ++i)
diff --git a/source/multigrid/mg_transfer_prebuilt.cc b/source/multigrid/mg_transfer_prebuilt.cc
index 84a3eb43d0..fa9d72f02e 100644
--- a/source/multigrid/mg_transfer_prebuilt.cc
+++ b/source/multigrid/mg_transfer_prebuilt.cc
@@ -150,6 +150,7 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
   // for a cell and one of its children
   std::vector<types::global_dof_index> dof_indices_parent (dofs_per_cell);
   std::vector<types::global_dof_index> dof_indices_child (dofs_per_cell);
+  std::vector<types::global_dof_index> entries (dofs_per_cell);
 
   // for each level: first build the sparsity
   // pattern of the matrices and then build the
@@ -158,20 +159,38 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
   // level which have children
   for (unsigned int level=0; level<n_levels-1; ++level)
     {
-
-      // reset the dimension of the structure.
-      // note that for the number of entries
-      // per row, the number of parent dofs
-      // coupling to a child dof is
-      // necessary. this, of course, is the
-      // number of degrees of freedom per
+      // find the locally relevant level dofs for setting the writable rows of
+      // the sparsity pattern
+      IndexSet level_p1_relevant_dofs = mg_dof.locally_owned_mg_dofs(level+1);
+      std::vector<types::global_dof_index> relevant_dofs;
+      for (typename DoFHandler<dim,spacedim>::cell_iterator cell=mg_dof.begin(level);
+           cell != mg_dof.end(level); ++cell)
+        if (cell->has_children() &&
+            ( mg_dof.get_tria().locally_owned_subdomain()==numbers::invalid_subdomain_id
+              || cell->level_subdomain_id()==mg_dof.get_tria().locally_owned_subdomain()
+            ))
+          for (unsigned int child=0; child<cell->n_children(); ++child)
+            {
+              cell->child(child)->get_mg_dof_indices (dof_indices_child);
+              for (unsigned int i=0; i<dofs_per_cell; ++i)
+                if (!level_p1_relevant_dofs.is_element(dof_indices_child[i]))
+                  relevant_dofs.push_back(dof_indices_child[i]);
+            }
+      std::sort(relevant_dofs.begin(), relevant_dofs.end());
+      level_p1_relevant_dofs.add_indices(relevant_dofs.begin(),
+                                         std::unique(relevant_dofs.begin(),
+                                                     relevant_dofs.end()));
+
+      // reset the dimension of the structure.  note that for the number of
+      // entries per row, the number of parent dofs coupling to a child dof is
+      // necessary. this, of course, is the number of degrees of freedom per
       // cell
-      // increment dofs_per_cell
-      // since a useless diagonal
-      // element will be stored
+      //
+      // increment dofs_per_cell since a useless diagonal element will be
+      // stored
       DynamicSparsityPattern dsp (sizes[level+1],
-                                  sizes[level]);
-      std::vector<types::global_dof_index> entries (dofs_per_cell);
+                                  sizes[level],
+                                  level_p1_relevant_dofs);
       for (typename DoFHandler<dim,spacedim>::cell_iterator cell=mg_dof.begin(level);
            cell != mg_dof.end(level); ++cell)
         if (cell->has_children() &&
@@ -316,13 +335,13 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
 
               if (global_mine && level_mine)
                 copy_indices[level].push_back(
-                  std::pair<unsigned int, unsigned int> (global_dof_indices[i], level_dof_indices[i]));
+                  std::make_pair (global_dof_indices[i], level_dof_indices[i]));
               else if (level_mine)
                 copy_indices_from_me[level].push_back(
-                  std::pair<unsigned int, unsigned int> (global_dof_indices[i], level_dof_indices[i]));
+                  std::make_pair (global_dof_indices[i], level_dof_indices[i]));
               else if (global_mine)
                 copy_indices_to_me[level].push_back(
-                  std::pair<unsigned int, unsigned int> (global_dof_indices[i], level_dof_indices[i]));
+                  std::make_pair (global_dof_indices[i], level_dof_indices[i]));
               else
                 continue;
 
@@ -334,7 +353,7 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
   // If we are in debugging mode, we order the copy indices, so we get
   // more reliable output for regression texts
 #ifdef DEBUG
-  std::less<std::pair<types::global_dof_index, unsigned int> > compare;
+  std::less<std::pair<types::global_dof_index, types::global_dof_index> > compare;
   for (unsigned int level=0; level<copy_indices.size(); ++level)
     std::sort(copy_indices[level].begin(), copy_indices[level].end(), compare);
   for (unsigned int level=0; level<copy_indices_from_me.size(); ++level)
diff --git a/source/multigrid/mg_transfer_prebuilt.inst.in b/source/multigrid/mg_transfer_prebuilt.inst.in
index fcf4a703f4..555809d410 100644
--- a/source/multigrid/mg_transfer_prebuilt.inst.in
+++ b/source/multigrid/mg_transfer_prebuilt.inst.in
@@ -28,7 +28,7 @@ for (deal_II_dimension : DIMENSIONS; V1 : VECTORS_WITH_MATRIX)
   }
 
 for (deal_II_dimension : DIMENSIONS; V1,V2 : DEAL_II_VEC_TEMPLATES; S1, S2 : REAL_SCALARS)
-  {    
+  {
     template void
       MGTransferPrebuilt<V1<S1> >::copy_to_mg (
 	const DoFHandler<deal_II_dimension>&, MGLevelObject<V1<S1> >&, const V2<S2>&) const;
@@ -40,8 +40,21 @@ for (deal_II_dimension : DIMENSIONS; V1,V2 : DEAL_II_VEC_TEMPLATES; S1, S2 : REA
 						 const MGLevelObject<V1<S1> >&) const;
   }
 
+for (deal_II_dimension : DIMENSIONS; S2 : REAL_SCALARS)
+  {
+    template void
+      MGTransferPrebuilt<parallel::distributed::Vector<double> >::copy_to_mg (
+	const DoFHandler<deal_II_dimension>&, MGLevelObject<parallel::distributed::Vector<double> >&, const parallel::distributed::Vector<S2>&) const;
+    template void
+      MGTransferPrebuilt<parallel::distributed::Vector<double> >::copy_from_mg (const DoFHandler<deal_II_dimension>&, parallel::distributed::Vector<S2>&,
+					     const MGLevelObject<parallel::distributed::Vector<double> >&) const;
+    template void
+      MGTransferPrebuilt<parallel::distributed::Vector<double> >::copy_from_mg_add (const DoFHandler<deal_II_dimension>&, parallel::distributed::Vector<S2>&,
+						 const MGLevelObject<parallel::distributed::Vector<double> >&) const;
+  }
+
 for(deal_II_dimension : DIMENSIONS)
-  {    
+  {
 #ifdef DEAL_II_WITH_TRILINOS
 
     template void
@@ -55,4 +68,3 @@ for(deal_II_dimension : DIMENSIONS)
 						 const MGLevelObject<TrilinosWrappers::MPI::Vector>&) const;
 #endif
   }
-