Fixed: FETools::inteprolation_difference was not working with TrilinosWrappers::MPI...

git-svn-id: https://svn.dealii.org/trunk@31686 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 4e97f5aaa3c85aa267325f104f1c9dfb5f562818..d699192944df975b125a1de9400d8389a0886255 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -211,6 +211,15 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: dealii::FETools::interpolation_difference was
+  not working for TrilinosWrappers::MPI::Vectors with ghost
+  entries. The TrilinosWrappers::VectorBase class has now a
+  get_mpi_communicator method similar to the PETSc vector
+  classes.
+  <br>
+  (Martin Steigemann, Martin Kronbichler, 2013/11/17)
+  </li>
+
   <li> Fixed: Bundled fparser is now compiled with FP_USE_THREAD_SAFE_EVAL in
   case of enabled threading support so that it is thread safe.
   <br>

diff --git a/deal.II/include/deal.II/lac/trilinos_vector_base.h b/deal.II/include/deal.II/lac/trilinos_vector_base.h
index a721cdce8c95ebd6135697994e81a8b4d8ff4da6..84d3f8291819b76bda4cb9e8dda3d3e43756cc51 100644 (file)
--- a/deal.II/include/deal.II/lac/trilinos_vector_base.h
+++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h
@@ -1051,6 +1051,13 @@ namespace TrilinosWrappers
      * consumption in bytes.
      */
     std::size_t memory_consumption () const;
+
+    /**
+     * Return a reference to the MPI
+     * communicator object in use with this
+     * object.
+     */
+    const MPI_Comm &get_mpi_communicator () const;
     //@}
 
     /**
@@ -2185,6 +2192,32 @@ namespace TrilinosWrappers
   }
 
 
+
+  inline
+  const MPI_Comm &
+  VectorBase::get_mpi_communicator () const
+  {
+    static MPI_Comm comm;
+
+#ifdef DEAL_II_WITH_MPI
+
+    const Epetra_MpiComm *mpi_comm
+      = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
+    comm = mpi_comm->Comm();
+
+#else
+
+    const Epetra_SerialComm *serial_comm
+      = dynamic_cast<const Epetra_SerialComm *>(&vector->Map().Comm());
+    comm = serial_comm->Comm();
+
+#endif
+
+    return comm;
+  }
+
+
+  
 #endif // DOXYGEN
 
 }

diff --git a/deal.II/source/fe/fe_tools_interpolate.cc b/deal.II/source/fe/fe_tools_interpolate.cc
index 8eef0c9e124eb6e95c4b62898729238aaec9dd1f..a1983255a1d9a746b081763c4b05fd415c5b35d6 100644 (file)
--- a/deal.II/source/fe/fe_tools_interpolate.cc
+++ b/deal.II/source/fe/fe_tools_interpolate.cc
@@ -93,15 +93,15 @@ namespace FETools
           // if u1 is a parallel distributed
           // PETSc vector, we check the local
           // size of u1 for safety
-          Assert(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size() == dof1.locally_owned_dofs().n_elements(),
-                 ExcDimensionMismatch(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size(), dof1.locally_owned_dofs().n_elements()));
+          Assert(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size() == dof1.n_locally_owned_dofs(),
+                 ExcDimensionMismatch(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size(), dof1.n_locally_owned_dofs()));
         };
 
     if (dynamic_cast<PETScWrappers::MPI::Vector *>(&u2) != 0)
       if (dynamic_cast<const DoFHandler<dim>*>(&dof2) != 0)
         {
-          Assert(dynamic_cast<PETScWrappers::MPI::Vector *>(&u2)->local_size() == dof2.locally_owned_dofs().n_elements(),
-                 ExcDimensionMismatch(dynamic_cast<PETScWrappers::MPI::Vector *>(&u2)->local_size(), dof2.locally_owned_dofs().n_elements()));
+          Assert(dynamic_cast<PETScWrappers::MPI::Vector *>(&u2)->local_size() == dof2.n_locally_owned_dofs(),
+                 ExcDimensionMismatch(dynamic_cast<PETScWrappers::MPI::Vector *>(&u2)->local_size(), dof2.n_locally_owned_dofs()));
         };
 #endif
     // allocate vectors at maximal
@@ -254,15 +254,15 @@ namespace FETools
           // if u1 is a parallel distributed
           // PETSc vector, we check the local
           // size of u1 for safety
-          Assert(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size() == dof1.locally_owned_dofs().n_elements(),
-                 ExcDimensionMismatch(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size(), dof1.locally_owned_dofs().n_elements()));
+          Assert(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size() == dof1.n_locally_owned_dofs(),
+                 ExcDimensionMismatch(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size(), dof1.n_locally_owned_dofs()));
         };
 
     if (dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated) != 0)
       if (dynamic_cast<const DoFHandler<dim>*>(&dof1) != 0)
         {
-          Assert(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated)->local_size() == dof1.locally_owned_dofs().n_elements(),
-                 ExcDimensionMismatch(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated)->local_size(), dof1.locally_owned_dofs().n_elements()));
+          Assert(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated)->local_size() == dof1.n_locally_owned_dofs(),
+                 ExcDimensionMismatch(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated)->local_size(), dof1.n_locally_owned_dofs()));
         };
 #endif
 
@@ -334,15 +334,15 @@ namespace FETools
           // if u1 is a parallel distributed
           // PETSc vector, we check the local
           // size of u1 for safety
-          Assert(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size() == dof1.locally_owned_dofs().n_elements(),
-                 ExcDimensionMismatch(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size(), dof1.locally_owned_dofs().n_elements()));
+          Assert(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size() == dof1.n_locally_owned_dofs(),
+                 ExcDimensionMismatch(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size(), dof1.n_locally_owned_dofs()));
         };
 
     if (dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated) != 0)
       if (dynamic_cast<const DoFHandler<dim>*>(&dof1) != 0)
         {
-          Assert(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated)->local_size() == dof1.locally_owned_dofs().n_elements(),
-                 ExcDimensionMismatch(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated)->local_size(), dof1.locally_owned_dofs().n_elements()));
+          Assert(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated)->local_size() == dof1.n_locally_owned_dofs(),
+                 ExcDimensionMismatch(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_interpolated)->local_size(), dof1.n_locally_owned_dofs()));
         };
 #endif
 
@@ -449,12 +449,41 @@ namespace FETools
         DoFTools::extract_locally_relevant_dofs (dof2,
                                                  dof2_locally_relevant_dofs);
 
-        PETScWrappers::MPI::Vector  u2_out (u1.get_mpi_communicator(),
-                                            dof2_locally_owned_dofs);
+        PETScWrappers::MPI::Vector  u2_out (dof2_locally_owned_dofs,
+                                            u1.get_mpi_communicator());
         interpolate(dof1, u1, dof2, constraints2, u2_out);
-        PETScWrappers::MPI::Vector  u2 (u1.get_mpi_communicator(),
-                                        dof2_locally_owned_dofs,
-                                        dof2_locally_relevant_dofs);
+        PETScWrappers::MPI::Vector  u2 (dof2_locally_owned_dofs,
+                                        dof2_locally_relevant_dofs,
+                                        u1.get_mpi_communicator());
+        u2 = u2_out;
+        interpolate(dof2, u2, dof1, constraints1, u1_interpolated);
+      }
+#endif
+
+      // special version for Trilinos
+#ifdef DEAL_II_USE_TRILINOS
+      template <int dim, int spacedim>
+      void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
+                             const ConstraintMatrix &constraints1,
+                             const TrilinosWrappers::MPI::Vector &u1,
+                             const DoFHandler<dim,spacedim> &dof2,
+                             const ConstraintMatrix &constraints2,
+                             TrilinosWrappers::MPI::Vector &u1_interpolated)
+      {
+        // if u1 is a parallel distributed Trilinos vector, we create a
+        // vector u2 with based on the sets of locally owned and relevant
+        // dofs of dof2
+        IndexSet  dof2_locally_owned_dofs = dof2.locally_owned_dofs();
+        IndexSet  dof2_locally_relevant_dofs;
+        DoFTools::extract_locally_relevant_dofs (dof2,
+                                                 dof2_locally_relevant_dofs);
+
+        TrilinosWrappers::MPI::Vector  u2_out (dof2_locally_owned_dofs,
+                                               u1.get_mpi_communicator());
+        interpolate(dof1, u1, dof2, constraints2, u2_out);
+        TrilinosWrappers::MPI::Vector  u2 (dof2_locally_owned_dofs,
+                                           dof2_locally_relevant_dofs,
+                                           u1.get_mpi_communicator());
         u2 = u2_out;
         interpolate(dof2, u2, dof1, constraints1, u1_interpolated);
       }
@@ -537,15 +566,15 @@ namespace FETools
           // if u1 is a parallel distributed
           // PETSc vector, we check the local
           // size of u1 for safety
-          Assert(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size() == dof1.locally_owned_dofs().n_elements(),
-                 ExcDimensionMismatch(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size(), dof1.locally_owned_dofs().n_elements()));
+          Assert(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size() == dof1.n_locally_owned_dofs(),
+                 ExcDimensionMismatch(dynamic_cast<const PETScWrappers::MPI::Vector *>(&u1)->local_size(), dof1.n_locally_owned_dofs()));
         };
 
     if (dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_difference) != 0)
       if (dynamic_cast<const DoFHandler<dim>*>(&dof1) != 0)
         {
-          Assert(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_difference)->local_size() == dof1.locally_owned_dofs().n_elements(),
-                 ExcDimensionMismatch(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_difference)->local_size(), dof1.locally_owned_dofs().n_elements()));
+          Assert(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_difference)->local_size() == dof1.n_locally_owned_dofs(),
+                 ExcDimensionMismatch(dynamic_cast<PETScWrappers::MPI::Vector *>(&u1_difference)->local_size(), dof1.n_locally_owned_dofs()));
         };
 #endif
 
@@ -598,6 +627,50 @@ namespace FETools
 
 
 
+  namespace internal
+  {
+    namespace
+    {
+      template <int dim, class InVector, class OutVector, int spacedim>
+      void interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
+                                     const ConstraintMatrix &constraints1,
+                                     const InVector &u1,
+                                     const DoFHandler<dim,spacedim> &dof2,
+                                     const ConstraintMatrix &constraints2,
+                                     OutVector &u1_difference)
+      {
+        back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference);
+        u1_difference.sadd(-1, u1);
+      }
+
+      // special version for Trilinos
+#ifdef DEAL_II_USE_TRILINOS
+      template <int dim, int spacedim>
+      void interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
+                                     const ConstraintMatrix &constraints1,
+                                     const TrilinosWrappers::MPI::Vector &u1,
+                                     const DoFHandler<dim,spacedim> &dof2,
+                                     const ConstraintMatrix &constraints2,
+                                     TrilinosWrappers::MPI::Vector &u1_difference)
+      {
+        back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference);
+
+        // Trilinos vectors with and without ghost entries are very different
+        // and we cannot use the sadd function directly, so we have to create
+        // a completely distributed vector first and copy the local entries
+        // from the vector with ghost entries
+        TrilinosWrappers::MPI::Vector  u1_completely_distributed (u1_difference.vector_partitioner ());
+
+        u1_completely_distributed = u1;
+
+        u1_difference.sadd(-1, u1_completely_distributed);
+      }
+#endif
+    }
+  }
+
+
+
   template <int dim, class InVector, class OutVector, int spacedim>
   void interpolation_difference(const DoFHandler<dim,spacedim> &dof1,
                                 const ConstraintMatrix &constraints1,
@@ -615,8 +688,7 @@ namespace FETools
       interpolation_difference(dof1, u1, dof2.get_fe(), u1_difference);
     else
       {
-        back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference);
-        u1_difference.sadd(-1, u1);
+        internal::interpolation_difference(dof1, constraints1, u1, dof2, constraints2, u1_difference);
       }
   }