]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use new features introduced in branch_distributed_grids even in mainline.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 1 Feb 2010 12:41:43 +0000 (12:41 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 1 Feb 2010 12:41:43 +0000 (12:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@20481 0785d39b-7218-0410-832d-ea1e28bc413d

13 files changed:
deal.II/base/include/base/index_set.h
deal.II/base/source/index_set.cc
deal.II/lac/include/lac/constraint_matrix.h
deal.II/lac/include/lac/constraint_matrix.templates.h
deal.II/lac/include/lac/petsc_solver.h
deal.II/lac/include/lac/trilinos_precondition.h
deal.II/lac/include/lac/trilinos_vector.h
deal.II/lac/source/constraint_matrix.cc
deal.II/lac/source/petsc_matrix_base.cc
deal.II/lac/source/petsc_solver.cc
deal.II/lac/source/trilinos_precondition.cc
deal.II/lac/source/trilinos_vector.cc
deal.II/lac/source/trilinos_vector_base.cc

index 3f4dc6b0bb303d42c7de1d2189be124141f2576c..2ca1086285a6e1e17d1707c7b82bef9f2859f816 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2009 by the deal.II authors
+//    Copyright (C) 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -150,6 +150,24 @@ class IndexSet
                                      */
     unsigned int index_within_set (const unsigned int n) const;
 
+                                    /**
+                                     * Each index set can be
+                                     * represented as the union of a
+                                     * number of contiguous intervals
+                                     * of indices, where if necessary
+                                     * intervals may only consist of
+                                     * individual elements to
+                                     * represent isolated members of
+                                     * the index set.
+                                     *
+                                     * This function returns the
+                                     * minimal number of such
+                                     * intervals that are needed to
+                                     * represent the index set under
+                                     * consideration.
+                                     */
+    unsigned int n_intervals () const;
+    
                                     /**
                                      * Compress the internal
                                      * representation by merging
@@ -206,6 +224,14 @@ class IndexSet
     IndexSet get_view (const unsigned int begin,
                       const unsigned int end) const;
 
+                                    /**
+                                     * Outputs a text representation of this
+                                     * IndexSet to the given stream. Used for
+                                     * testing.
+                                     */
+    template <class STREAM>
+    void print(STREAM &out) const;
+    
 #ifdef DEAL_II_USE_TRILINOS
                                     /**
                                      * Given an MPI communicator,
@@ -655,6 +681,26 @@ IndexSet::operator != (const IndexSet &is) const
   return ranges != is.ranges;
 }
 
+template <class STREAM>
+inline
+void
+IndexSet::print (STREAM &out) const
+{
+  compress();
+  out << "{";
+  std::vector<Range>::const_iterator p;
+  for (p = ranges.begin(); p != ranges.end(); ++p)
+    {
+      if (p->end-p->begin==1)
+       out << p->begin;
+      else
+       out << "[" << p->begin << "," << p->end-1 << "]";
+
+      if (p !=--ranges.end())
+       out << ", ";
+    }
+  out << "}" << std::endl;
+}
 
 
 DEAL_II_NAMESPACE_CLOSE
index 385571375b5d52738efffd2ecd5ecc40856251b9..ac5e5cf71f122926dfbc17684b9ad2fe0c406908 100644 (file)
@@ -154,6 +154,15 @@ IndexSet::operator & (const IndexSet &is) const
 
 
 
+unsigned int
+IndexSet::n_intervals () const
+{
+  compress ();
+  return ranges.size();
+}
+
+
+
 IndexSet
 IndexSet::get_view (const unsigned int begin,
                    const unsigned int end) const
index 789519f1567462f68dd9c4dd231e6d6f3f029ee1..4b7f379544c382b596ff8b21748c146bce22ef10 100644 (file)
@@ -1546,6 +1546,17 @@ class ConstraintMatrix : public Subscriptor
                    << "to another DoF with number " << arg1
                    << ", which however is constrained by this object. This is not"
                    << " allowed.");
+                                    /**
+                                     * Exception
+                                     *
+                                     * @ingroup Exceptions
+                                     */
+    DeclException1 (ExcRowNotStoredHere,
+                   int,
+                   << "The index set given to this constraint matrix indicates "
+                   << "constraints for degree of freedom " << arg1
+                   << " should not be stored by this object, but a constraint "
+                   << "is being added.");
 
   private:
 
@@ -2075,7 +2086,8 @@ ConstraintMatrix::calculate_line_index (const unsigned int line) const
   if (!local_lines.size())
     return line;
 
-  Assert(local_lines.is_element(line), ExcInternalError());
+  Assert(local_lines.is_element(line),
+        ExcRowNotStoredHere(line));
 
   return local_lines.index_within_set(line);
 }
index e956c626b43fc2b799b415bffe5b345e053dd9ad..37aa2128da7e5134d930871d0b3db46cb303cd58 100644 (file)
@@ -1077,7 +1077,7 @@ namespace internals
                                    // specialized sort and insert functions.
   struct GlobalRowsFromLocal
   {
-    GlobalRowsFromLocal (const unsigned int n_local_rows) 
+    GlobalRowsFromLocal (const unsigned int n_local_rows)
       :
       total_row_indices (n_local_rows), n_active_rows (n_local_rows),
       n_inhomogeneous_rows (0){};
@@ -1110,9 +1110,9 @@ namespace internals
                [index_in_constraint]).second; };
     bool have_indirect_rows () const { return data_cache.element_size; };
     void insert_constraint (const unsigned int constrained_local_dof)
-      { --n_active_rows; 
+      { --n_active_rows;
        total_row_indices[n_active_rows].local_row = constrained_local_dof; }
-    unsigned int last_constrained_local_row () 
+    unsigned int last_constrained_local_row ()
       { Assert (total_row_indices.back().global_row == numbers::invalid_unsigned_int,
                ExcInternalError());
        return total_row_indices.back().local_row; };
@@ -1447,7 +1447,7 @@ namespace internals
   {
     if (value != 0.)
       {
-       Assert (col_ptr != 0, 
+       Assert (col_ptr != 0,
                typename SparseMatrix<number>::ExcInvalidIndex (row, column));
        if (are_on_diagonal)
          {
@@ -1491,7 +1491,7 @@ namespace internals
     const unsigned int loc_row = global_rows.local_row(i);
 
 #ifdef DEBUG
-    const unsigned int * col_ptr = sparsity.n_nonzero_elements() == 0 ? 0 : 
+    const unsigned int * col_ptr = sparsity.n_nonzero_elements() == 0 ? 0 :
       &sparsity_struct[row_start[row]];
     number * val_ptr = sparsity.n_nonzero_elements() == 0 ? 0 :
       &sparse_matrix->global_entry (row_start[row]);
index 0aac5ad434fa590d6eb9fa97648bb14e75f64b58..545ce0d0ba40dcc60f68bdd8ca839f7e6594e556 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007, 2009 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -42,6 +42,12 @@ namespace PETScWrappers
  * classes simply set the right flags to select one solver or another, or to
  * set certain parameters for individual solvers.
  *
+ * Note: Repeated calls to solve() on a solver object with a Preconditioner
+ * must be used with care. The preconditioner is initialized in the first call
+ * to solve() and subsequent calls reuse the solver and preconditioner
+ * object. This is done for performance reasons. The solver and preconditioner
+ * can be reset by calling reset().
+ *
  * One of the gotchas of PETSc is that -- in particular in MPI mode -- it
  * often does not produce very helpful error messages. In order to save
  * other users some time in searching a hard to track down error, here is
@@ -101,7 +107,10 @@ namespace PETScWrappers
                                         * classes and the object passed as a
                                         * preconditioner, one of the linear
                                         * solvers and preconditioners of PETSc
-                                        * is chosen.
+                                        * is chosen.  Repeated calls to
+                                        * solve() do not reconstruct the
+                                        * preconditioner for performance
+                                        * reasons. See class Documentation.
                                         */
       void
       solve (const MatrixBase         &A,
@@ -110,6 +119,14 @@ namespace PETScWrappers
              const PreconditionerBase &preconditioner);
 
 
+                                      /**
+                                       * Resets the contained preconditioner
+                                       * and solver object. See class
+                                       * description for more details.
+                                       */
+      virtual void reset();
+                                       
+      
                                        /**
                                         * Access to object that controls
                                         * convergence.
index edfd3415cdfb4f8b9406949d941cfed54517761f..8f25060460386da1a3e58607eed8e72b85ce6cb5 100644 (file)
@@ -1308,9 +1308,11 @@ namespace TrilinosWrappers
                                        */
        AdditionalData (const bool                             elliptic = true,
                        const bool                             higher_order_elements = false,
+                       const unsigned int                     n_cycles = 1,
+                       const bool                             w_cyle = false,
                        const double                           aggregation_threshold = 1e-4,
                        const std::vector<std::vector<bool> > &constant_modes = std::vector<std::vector<bool> > (1),
-                       const unsigned int                     smoother_sweeps = 3,
+                       const unsigned int                     smoother_sweeps = 2,
                        const unsigned int                     smoother_overlap = 0,
                        const bool                             output_details = false);
 
@@ -1335,6 +1337,20 @@ namespace TrilinosWrappers
                                        */
        bool higher_order_elements;
 
+                                      /**
+                                       * Defines how many multigrid cycles
+                                       * should be performed by the
+                                       * preconditioner.
+                                       */
+       unsigned int n_cycles;
+
+                                      /**
+                                       * Defines whether a w-cycle should be
+                                       * used instead of the standard setting
+                                       * of a v-cycle.
+                                       */
+       bool w_cycle;
+
                                       /**
                                        * This threshold tells the AMG setup
                                        * how the coarsening should be
index ec0a4e82e79602ac46e82ed65757192bf5ec75b4..ae24942466b7f42cd7560f7f80df667073131105 100644 (file)
@@ -51,6 +51,7 @@ namespace TrilinosWrappers
  */
   namespace MPI
   {
+    class BlockVector;
 
 /**
  * This class implements a wrapper to use the Trilinos distributed
@@ -232,6 +233,9 @@ namespace TrilinosWrappers
                     const bool        fast = false,
                     const bool        allow_different_maps = false);
 
+       void reinit (const BlockVector &v,
+                    const bool         import_data = false);
+
                                        /**
                                        * Set all components of the
                                        * vector to the given number @p
index 4969878ad6dbcea57bcb3e5b697b9f8fd6c87b56..a4da438d13084497beca0d021c1d41f2a545fd26 100644 (file)
@@ -264,10 +264,19 @@ void ConstraintMatrix::close ()
                                           // this line (including
                                           // ones that we have
                                           // appended in this go
-                                          // around)
+                                          // around) and see whether
+                                          // they are further
+                                          // constrained. ignore
+                                          // elements that we don't
+                                          // store on the current
+                                          // processor
          unsigned int entry = 0;
          while (entry < line->entries.size())
-           if (is_constrained (line->entries[entry].first))
+           if (((local_lines.size() == 0)
+                ||
+                (local_lines.is_element(line->entries[entry].first)))
+               &&
+               is_constrained (line->entries[entry].first))
              {
                                                 // ok, this entry is
                                                 // further
@@ -504,20 +513,27 @@ void ConstraintMatrix::close ()
 #ifdef DEBUG
                                   // if in debug mode: check that no
                                   // dof is constraint to another dof
-                                  // that is also constrained
+                                  // that is also
+                                  // constrained. exclude dofs from
+                                  // this check whose constraint
+                                  // lines would not be store on the
+                                  // local processor
   for (std::vector<ConstraintLine>::const_iterator line=lines.begin();
        line!=lines.end(); ++line)
     for (std::vector<std::pair<unsigned int,double> >::const_iterator
           entry=line->entries.begin();
         entry!=line->entries.end(); ++entry)
-      {
-                                        // make sure that
-                                        // entry->first is not the
-                                        // index of a line itself
-       const bool is_circle = is_constrained(entry->first);
-       Assert (is_circle == false,
-               ExcDoFConstrainedToConstrainedDoF(line->line, entry->first));
-      };
+      if ((local_lines.size() == 0)
+         ||
+         (local_lines.is_element(entry->first)))
+       {
+                                          // make sure that
+                                          // entry->first is not the
+                                          // index of a line itself
+         const bool is_circle = is_constrained(entry->first);
+         Assert (is_circle == false,
+                 ExcDoFConstrainedToConstrainedDoF(line->line, entry->first));
+       }
 #endif
 
   sorted = true;
@@ -1919,6 +1935,7 @@ void ConstraintMatrix::condense (BlockCompressedSimpleSparsityPattern &sparsity)
                                   // on the local range of the vector. But
                                   // we need access to values that the
                                   // local nodes are constrained to.
+
 template<>
 void
 ConstraintMatrix::distribute (TrilinosWrappers::MPI::Vector &vec) const
@@ -1992,6 +2009,98 @@ ConstraintMatrix::distribute (TrilinosWrappers::MPI::Vector &vec) const
     }
 }
 
+
+
+template<>
+void
+ConstraintMatrix::distribute (TrilinosWrappers::MPI::BlockVector &vec) const
+{
+  IndexSet my_indices (vec.size());
+  for (unsigned int block=0; block<vec.n_blocks(); ++block)
+    {
+      typedef std::vector<ConstraintLine>::const_iterator constraint_iterator;
+      ConstraintLine index_comparison;
+      index_comparison.line = vec.block(block).local_range().first
+       +vec.get_block_indices().block_start(block);
+      const constraint_iterator begin_my_constraints =
+       std::lower_bound (lines.begin(),lines.end(),index_comparison);
+
+      index_comparison.line = vec.block(block).local_range().second
+       +vec.get_block_indices().block_start(block);
+
+      const constraint_iterator end_my_constraints
+       = std::lower_bound(lines.begin(),lines.end(),index_comparison);
+
+                                  // Here we search all the indices that we
+                                  // need to have read-access to - the local
+                                  // nodes and all the nodes that the
+                                  // constraints indicate. No caching done
+                                  // yet. would need some more clever data
+                                  // structures for doing that.
+      const std::pair<unsigned int, unsigned int>
+       local_range = vec.block(block).local_range();
+
+      my_indices.add_range (local_range.first, local_range.second);
+
+      std::set<unsigned int> individual_indices;
+      for (constraint_iterator it = begin_my_constraints;
+          it != end_my_constraints; ++it)
+       for (unsigned int i=0; i<it->entries.size(); ++i)
+         if ((it->entries[i].first < local_range.first)
+             ||
+             (it->entries[i].first >= local_range.second))
+           individual_indices.insert (it->entries[i].first);
+
+      my_indices.add_indices (individual_indices.begin(),
+                             individual_indices.end());
+    }
+
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+      const Epetra_MpiComm *mpi_comm
+       = dynamic_cast<const Epetra_MpiComm*>(&vec.block(0).trilinos_vector().Comm());
+
+      Assert (mpi_comm != 0, ExcInternalError());
+
+      TrilinosWrappers::MPI::Vector vec_distribute 
+       (my_indices.make_trilinos_map (mpi_comm->Comm(), true));
+#else
+      TrilinosWrappers::MPI::Vector vec_distribute 
+       (my_indices.make_trilinos_map (MPI_COMM_WORLD, true)));
+#endif
+
+                                  // here we import the data
+  vec_distribute.reinit(vec,true);
+
+  for (unsigned int block=0; block<vec.n_blocks(); ++block)
+    {
+      typedef std::vector<ConstraintLine>::const_iterator constraint_iterator;
+      ConstraintLine index_comparison;
+      index_comparison.line = vec.block(block).local_range().first
+       +vec.get_block_indices().block_start(block);
+      const constraint_iterator begin_my_constraints =
+       std::lower_bound (lines.begin(),lines.end(),index_comparison);
+
+      index_comparison.line = vec.block(block).local_range().second
+       +vec.get_block_indices().block_start(block);
+
+      const constraint_iterator end_my_constraints
+       = std::lower_bound(lines.begin(),lines.end(),index_comparison);
+
+      for (constraint_iterator it = begin_my_constraints;
+          it != end_my_constraints; ++it)
+       {
+                                      // fill entry in line
+                                      // next_constraint.line by adding the
+                                      // different contributions
+         double new_value = it->inhomogeneity;
+         for (unsigned int i=0; i<it->entries.size(); ++i)
+           new_value += (vec_distribute(it->entries[i].first) *
+                         it->entries[i].second);
+         vec(it->line) = new_value;
+       }
+    }
+}
+
 #endif
 
 
index 868e7261b012631928c4f9aaf0659df8a7953145..f19b3944e800b289104180c7a0b9a24aa894c78b 100644 (file)
@@ -620,7 +620,7 @@ namespace PETScWrappers
        std::cout << "*** PETSC-Matrix: num-allocs = "
                  << info.mallocs << " ***" << std::endl;
 
-      return sizeof(*this) + info.memory;
+      return sizeof(*this) + static_cast<unsigned int>(info.memory);
   }
 
 }
index cc9bf46911cbd601c975371e683cbf46c62919b7..94ec618415787ad664bf37286ec0d09961af5f84 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2006, 2008, 2009 by the deal.II authors
+//    Copyright (C) 2004, 2006, 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -198,8 +198,8 @@ namespace PETScWrappers
 #  endif
 
 #endif
-                                     // destroy solver object
-    solver_data.reset ();
+                                     // do not destroy solver object
+//    solver_data.reset ();
 
                                      // in case of failure: throw
                                      // exception
@@ -211,6 +211,14 @@ namespace PETScWrappers
 
 
 
+  void
+  SolverBase::reset()
+  {
+    solver_data.reset ();
+  }
+
+  
+
   SolverControl &
   SolverBase::control() const
   {
index f0f337bbff423da8c247190c7df8c7cff7ff04b1..2be68d6be1267b6898f1794c51013bed4e8b8e74 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2008 by the deal.II authors
+//    Copyright (C) 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be distributed
 //    without copyright and license information. Please refer
@@ -488,6 +488,8 @@ namespace TrilinosWrappers
   PreconditionAMG::AdditionalData::
   AdditionalData (const bool                             elliptic,
                  const bool                             higher_order_elements,
+                 const unsigned int                     n_cycles,
+                 const bool                             w_cycle,
                  const double                           aggregation_threshold,
                  const std::vector<std::vector<bool> > &constant_modes,
                  const unsigned int                     smoother_sweeps,
@@ -496,6 +498,8 @@ namespace TrilinosWrappers
                   :
                   elliptic (elliptic),
                  higher_order_elements (higher_order_elements),
+                 n_cycles (n_cycles),
+                 w_cycle (w_cycle),
                  aggregation_threshold (aggregation_threshold),
                  constant_modes (constant_modes),
                  smoother_sweeps (smoother_sweeps),
@@ -513,8 +517,6 @@ namespace TrilinosWrappers
     multilevel_operator.release();
 
     const unsigned int n_rows = matrix.m();
-    const unsigned int constant_modes_dimension =
-      additional_data.constant_modes.size();
 
                                        // Build the AMG preconditioner.
     Teuchos::ParameterList parameter_list;
@@ -559,6 +561,13 @@ namespace TrilinosWrappers
 
     parameter_list.set("smoother: sweeps",
                       static_cast<int>(additional_data.smoother_sweeps));
+    parameter_list.set("cycle applications", 
+                      static_cast<int>(additional_data.n_cycles));
+    if (additional_data.w_cycle == true)
+      parameter_list.set("prec type", "MGW");
+    else
+      parameter_list.set("prec type", "MGV");
+      
     parameter_list.set("smoother: Chebyshev alpha",10.);
     parameter_list.set("smoother: ifpack overlap",
                       static_cast<int>(additional_data.smoother_overlap));
@@ -572,23 +581,26 @@ namespace TrilinosWrappers
 
     const Epetra_Map & domain_map = matrix.domain_partitioner();
 
+    const unsigned int constant_modes_dimension =
+      additional_data.constant_modes.size();
     Epetra_MultiVector distributed_constant_modes (domain_map,
                                                   constant_modes_dimension);
 
     if (constant_modes_dimension > 1)
       {
-       Assert (n_rows == additional_data.constant_modes[0].size(),
-               ExcDimensionMismatch(n_rows,
-                                    additional_data.constant_modes[0].size()));
+       const bool constant_modes_are_global =
+         additional_data.constant_modes[0].size() == n_rows;
+       const unsigned int n_relevant_rows =
+         constant_modes_are_global ? n_rows : additional_data.constant_modes[0].size();
+       const unsigned int my_size = domain_map.NumMyElements();
+       if (constant_modes_are_global == false)
+         Assert (n_relevant_rows == my_size,
+                 ExcDimensionMismatch(n_relevant_rows, my_size));
        Assert (n_rows ==
                static_cast<unsigned int>(distributed_constant_modes.GlobalLength()),
                ExcDimensionMismatch(n_rows,
                                     distributed_constant_modes.GlobalLength()));
 
-       const unsigned int my_size = domain_map.NumMyElements();
-       Assert (my_size == static_cast<unsigned int>(domain_map.MaxLID()+1),
-               ExcDimensionMismatch (my_size, domain_map.MaxLID()+1));
-
                                        // Reshape null space as a
                                        // contiguous vector of
                                        // doubles so that Trilinos
@@ -596,7 +608,7 @@ namespace TrilinosWrappers
        for (unsigned int d=0; d<constant_modes_dimension; ++d)
          for (unsigned int row=0; row<my_size; ++row)
            {
-             int global_row_id = domain_map.GID(row);
+             int global_row_id = constant_modes_are_global ? domain_map.GID(row) : row;
              distributed_constant_modes.ReplaceMyValue(row, d,
                static_cast<double>(additional_data.constant_modes[d][global_row_id]));
            }
index 00795b8c23ff7545957df7e57100af039dc66cc2..ed08c205eaf585251b7bd4b47aebf0323552c5ce 100644 (file)
 #ifdef DEAL_II_USE_TRILINOS
 
 #  include <lac/trilinos_sparse_matrix.h>
+#  include <lac/trilinos_block_vector.h>
 #  include <cmath>
 #  include <Epetra_Import.h>
+#  include <Epetra_Vector.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -209,6 +211,73 @@ namespace TrilinosWrappers
 
 
 
+    void
+    Vector::reinit (const BlockVector &v,
+                   const bool         import_data)
+    {
+                                       // In case we do not allow to
+                                       // have different maps, this
+                                       // call means that we have to
+                                       // reset the vector. So clear
+                                       // the vector, initialize our
+                                       // map with the map in v, and
+                                       // generate the vector.
+      if (v.n_blocks() == 0)
+       return;
+
+                               // create a vector that holds all the elements
+                               // contained in the block vector. need to
+                               // manually create an Epetra_Map.
+      unsigned int n_elements = 0, added_elements = 0, block_offset = 0;
+      for (unsigned int block=0; block<v.n_blocks();++block)
+       n_elements += v.block(block).local_size();
+      std::vector<int> global_ids (n_elements, -1);
+      for (unsigned int block=0; block<v.n_blocks();++block)
+       {
+         int * glob_elements = v.block(block).vector_partitioner().MyGlobalElements();
+         for (unsigned int i=0; i<v.block(block).local_size(); ++i)
+           global_ids[added_elements++] = glob_elements[i] + block_offset;
+         block_offset += v.block(block).size();
+       }
+
+      Assert (n_elements == added_elements, ExcInternalError());
+      Epetra_Map new_map (v.size(), n_elements, &global_ids[0], 0,
+                         v.block(0).vector_partitioner().Comm());
+
+      if (import_data == false)
+       vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector (new_map));
+
+      Teuchos::RCP<Epetra_FEVector> actual_vec = (import_data == true) ? 
+       Teuchos::rcp (new Epetra_FEVector (new_map), true) : 
+       Teuchos::rcp (&*vector, false);
+
+      TrilinosScalar* entries = (*actual_vec)[0];
+      block_offset = 0;
+      for (unsigned int block=0; block<v.n_blocks();++block)
+       {
+         v.block(block).trilinos_vector().ExtractCopy (entries, 0);
+         entries += v.block(block).local_size();
+       }
+
+      if (import_data == true)
+        {
+         AssertThrow (static_cast<unsigned int>(actual_vec->GlobalLength()) 
+                      == v.size(),
+                      ExcDimensionMismatch (actual_vec->GlobalLength(), 
+                                            v.size()));
+
+         Epetra_Import data_exchange (vector->Map(), actual_vec->Map());
+
+         const int ierr = vector->Import(*actual_vec, data_exchange, Insert);
+         AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+         last_action = Insert;
+       }
+
+    }
+
+
+
     Vector &
     Vector::operator = (const Vector &v)
     {
@@ -226,8 +295,9 @@ namespace TrilinosWrappers
        {
          Assert (vector->Map().SameAs(v.vector->Map()) == true,
                  ExcMessage ("The Epetra maps in the assignment operator ="
-                             " do not match, even though the local_range "
-                             " seems to be the same. Check vector setup!"));
+                             " do not match, even though the local_range"
+                             " seems to be the same. Check vector setup or"
+                             " use reinit()!"));
 
          const int ierr = vector->Update(1.0, *v.vector, 0.0);
          Assert (ierr == 0, ExcTrilinosError(ierr));
@@ -240,7 +310,21 @@ namespace TrilinosWrappers
        {
          Epetra_Import data_exchange (vector->Map(), v.vector->Map());
          const int ierr = vector->Import(*v.vector, data_exchange, Insert);
-         AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+         Assert (ierr == 0, ExcTrilinosError(ierr));
+       }
+      else if (size() == v.size() && local_range().first>=v.local_range().first &&
+              local_range().second<=v.local_range().second &&
+              vector->Map().UniqueGIDs())
+       {
+         for (unsigned int i=0; i<local_size(); ++i)
+           {
+             Assert (v.trilinos_vector().Map().MyGID(vector->Map().GID(i)),
+                     ExcMessage ("The right hand side vector does not contain "
+                                 "all local elements present in the left hand "
+                                 "side vector in the assignment operation =, "
+                                 "which is not allowed. Use reinit()."));
+             (*vector)[0][i] = v(vector->Map().GID(i));
+           }
        }
       else
        {
index 5fba85de9376328586abcaa6874cd1d7705c8153..2cfee040e69794ca3ec83adb2549ec2b117783e3 100644 (file)
@@ -272,17 +272,9 @@ namespace TrilinosWrappers
   {
                                      // get a representation of the vector and
                                      // loop over all the elements
-    TrilinosScalar *start_ptr;
-    int leading_dimension;
-    int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-                                      // TODO: This
-                                      // won't work in parallel like
-                                      // this. Find out a better way to
-                                      // this in that case.
+    TrilinosScalar *start_ptr = (*vector)[0];
     const TrilinosScalar *ptr  = start_ptr,
-                         *eptr = start_ptr + size();
+                         *eptr = start_ptr + local_size();
     bool flag = true;
     while (ptr != eptr)
       {

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.