Make ConstraintMatrix::distribute also work for block vectors. Needs a couple extra...

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

diff --git a/deal.II/include/deal.II/lac/constraint_matrix.templates.h b/deal.II/include/deal.II/lac/constraint_matrix.templates.h
index a1c17c953826df619dc3bf571d982991a37b3ff4..5758495937ecac4fab9dbb6a043ed161ad6e74c3 100644 (file)
--- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h
+++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h
@@ -988,69 +988,52 @@ namespace internal
 #ifdef DEAL_II_WITH_TRILINOS
     void
     import_vector_with_ghost_elements (const TrilinosWrappers::MPI::Vector &vec,
-                                       const IndexSet       &/*locally_owned_elements*/,
-                                       const IndexSet       &needed_elements,
-                                       TrilinosWrappers::MPI::Vector       &output)
+                                       const IndexSet                      &/*locally_owned_elements*/,
+                                       const IndexSet                      &needed_elements,
+                                       TrilinosWrappers::MPI::Vector       &output,
+                                       const internal::bool2type<false>     /*is_block_vector*/)
     {
       output.reinit (needed_elements,
                      dynamic_cast<const Epetra_MpiComm *>(&vec.vector_partitioner().Comm())->GetMpiComm());
       output = vec;
     }
-
-
-    void
-    import_vector_with_ghost_elements (const TrilinosWrappers::Vector &vec,
-                                       const IndexSet       &,
-                                       const IndexSet       &,
-                                       TrilinosWrappers::Vector       &output)
-    {
-      output.reinit (vec.size());
-      output = vec;
-    }
 #endif
 
 #ifdef DEAL_II_WITH_PETSC
     void
     import_vector_with_ghost_elements (const PETScWrappers::MPI::Vector &vec,
-                                       const IndexSet       &locally_owned_elements,
-                                       const IndexSet       &needed_elements,
-                                       PETScWrappers::MPI::Vector       &output)
+                                       const IndexSet                   &locally_owned_elements,
+                                       const IndexSet                   &needed_elements,
+                                       PETScWrappers::MPI::Vector       &output,
+                                       const internal::bool2type<false>  /*is_block_vector*/)
     {
       output.reinit (vec.get_mpi_communicator(), locally_owned_elements, needed_elements);
       output = vec;
     }
-
-
-
-    void
-    import_vector_with_ghost_elements (const PETScWrappers::Vector &vec,
-                                       const IndexSet       &,
-                                       const IndexSet       &,
-                                       PETScWrappers::Vector       &output)
-    {
-      output.reinit (vec.size());
-      output = vec;
-    }
 #endif
 
     template <typename number>
     void
     import_vector_with_ghost_elements (const parallel::distributed::Vector<number> &vec,
-                                       const IndexSet       &locally_owned_elements,
-                                       const IndexSet       &needed_elements,
-                                       parallel::distributed::Vector<number>       &output)
+                                       const IndexSet                              &locally_owned_elements,
+                                       const IndexSet                              &needed_elements,
+                                       parallel::distributed::Vector<number>       &output,
+                                       const internal::bool2type<false>             /*is_block_vector*/)
     {
       output.reinit (locally_owned_elements, needed_elements, vec.get_mpi_communicator());
       output = vec;
     }
 
 
-    template <typename number>
+    // all other vector non-block vector types are sequential and we should
+    // not have this function called at all -- so throw an exception
+    template <typename Vector>
     void
-    import_vector_with_ghost_elements (const dealii::Vector<number> &vec,
-                                       const IndexSet   &locally_owned_elements,
-                                       const IndexSet   &needed_elements,
-                                       dealii::Vector<number>       &output)
+    import_vector_with_ghost_elements (const Vector                     &/*vec*/,
+                                       const IndexSet                   &/*locally_owned_elements*/,
+                                       const IndexSet                   &/*needed_elements*/,
+                                       Vector                           &/*output*/,
+                                       const internal::bool2type<false>  /*is_block_vector*/)
     {
       Assert (false, ExcMessage ("We shouldn't even get here!"));
     }
@@ -1059,20 +1042,26 @@ namespace internal
     // for block vectors, simply dispatch to the individual blocks
     template <class VectorType>
     void
-    import_vector_with_ghost_elements (const BlockVectorBase<VectorType> &vec,
-                                       const IndexSet   &locally_owned_elements,
-                                       const IndexSet   &needed_elements,
-                                       BlockVectorBase<VectorType>       &output)
+    import_vector_with_ghost_elements (const VectorType                &vec,
+                                       const IndexSet                  &locally_owned_elements,
+                                       const IndexSet                  &needed_elements,
+                                       VectorType                      &output,
+                                       const internal::bool2type<true>  /*is_block_vector*/)
     {
+      output.reinit (vec.n_blocks());
+
       types::global_dof_index block_start = 0;
       for (unsigned int b=0; b<vec.n_blocks(); ++b)
         {
           import_vector_with_ghost_elements (vec.block(b),
                                              locally_owned_elements.get_view (block_start, block_start+vec.block(b).size()),
                                              needed_elements.get_view (block_start, block_start+vec.block(b).size()),
-                                             output.block(b));
+                                             output.block(b),
+                                             internal::bool2type<false>());
           block_start += vec.block(b).size();
         }
+
+      output.collect_sizes ();
     }
   }
 }
@@ -1131,7 +1120,8 @@ ConstraintMatrix::distribute (VectorType &vec) const
       VectorType ghosted_vector;
       internal::import_vector_with_ghost_elements (vec,
                                                    vec_owned_elements, needed_elements,
-                                                   ghosted_vector);
+                                                   ghosted_vector,
+                                                   internal::bool2type<IsBlockVector<VectorType>::value>());
 
       for (constraint_iterator it = lines.begin();
           it != lines.end(); ++it)

diff --git a/deal.II/include/deal.II/lac/petsc_block_vector.h b/deal.II/include/deal.II/lac/petsc_block_vector.h
index 788d8872b29a8f0899ca73d865f6d7ba9921418c..0c55816a93632833e9672a49ad83e48fe9188b8a 100644 (file)
--- a/deal.II/include/deal.II/lac/petsc_block_vector.h
+++ b/deal.II/include/deal.II/lac/petsc_block_vector.h
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007, 2010, 2012 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2010, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -259,6 +259,19 @@ namespace PETScWrappers
     void reinit (const BlockVector &V,
                  const bool         fast=false);
 
+    /**
+     * Change the number of blocks to
+     * <tt>num_blocks</tt>. The individual
+     * blocks will get initialized with
+     * zero size, so it is assumed that
+     * the user resizes the
+     * individual blocks by herself
+     * in an appropriate way, and
+     * calls <tt>collect_sizes</tt>
+     * afterwards.
+     */
+    void reinit (const unsigned int num_blocks);
+
     /**
      * Swap the contents of this
      * vector and the other vector
@@ -458,6 +471,15 @@ namespace PETScWrappers
 
 
 
+  inline
+  void
+  BlockVector::reinit (const unsigned int num_blocks)
+  {
+    reinit (num_blocks, 0, true);
+  }
+
+
+
   inline
   void
   BlockVector::swap (BlockVector &v)

diff --git a/deal.II/include/deal.II/lac/petsc_parallel_block_vector.h b/deal.II/include/deal.II/lac/petsc_parallel_block_vector.h
index 33482381d21fe4c46b3c11cf964b765ac061de9a..70c3606c64ca84e6b8d0a9e087ff183999727cab 100644 (file)
--- a/deal.II/include/deal.II/lac/petsc_parallel_block_vector.h
+++ b/deal.II/include/deal.II/lac/petsc_parallel_block_vector.h
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007, 2009, 2010, 2012 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2009, 2010, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -279,6 +279,19 @@ namespace PETScWrappers
       void reinit (const std::vector<IndexSet> &parallel_partitioning,
                    const MPI_Comm              &communicator);
 
+      /**
+       * Change the number of blocks to
+       * <tt>num_blocks</tt>. The individual
+       * blocks will get initialized with
+       * zero size, so it is assumed that
+       * the user resizes the
+       * individual blocks by herself
+       * in an appropriate way, and
+       * calls <tt>collect_sizes</tt>
+       * afterwards.
+       */
+      void reinit (const unsigned int num_blocks);
+
       /**
        * Return a reference to the MPI
        * communicator object in use with

diff --git a/deal.II/source/lac/petsc_parallel_block_vector.cc b/deal.II/source/lac/petsc_parallel_block_vector.cc
index fc11d85a10e26022c0cb0a568a315577da83abcb..49318b38400c5560da0be13af434c45210d046e1 100644 (file)
--- a/deal.II/source/lac/petsc_parallel_block_vector.cc
+++ b/deal.II/source/lac/petsc_parallel_block_vector.cc
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2008, 2012 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2008, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -35,6 +35,21 @@ namespace PETScWrappers
 
       return *this;
     }
+
+
+    void
+    BlockVector::reinit (const unsigned int num_blocks)
+    {
+      std::vector<unsigned int> block_sizes (num_blocks, 0);
+      this->block_indices.reinit (block_sizes);
+      if (this->components.size() != this->n_blocks())
+        this->components.resize(this->n_blocks());
+
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        components[i].reinit (MPI_Comm(), 0, 0);
+
+      collect_sizes();
+    }
   }
 
 }