]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merged from trunk.
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 5 Jan 2014 19:09:16 +0000 (19:09 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 5 Jan 2014 19:09:16 +0000 (19:09 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32170 0785d39b-7218-0410-832d-ea1e28bc413d

36 files changed:
deal.II/bundled/CMakeLists.txt
deal.II/bundled/boost-1.49.0/libs/iostreams/src/CMakeLists.txt
deal.II/cmake/modules/FindDEALII_LAPACK.cmake
deal.II/cmake/setup_cached_variables.cmake
deal.II/doc/doxygen/CMakeLists.txt
deal.II/doc/doxygen/headers/vector_valued.h
deal.II/doc/doxygen/options.dox.in
deal.II/doc/news/changes.h
deal.II/doc/publications/dealii.bib
deal.II/doc/publications/index.html
deal.II/examples/step-12/step-12.cc
deal.II/include/deal.II/lac/block_sparsity_pattern.h
deal.II/include/deal.II/lac/relaxation_block.h
deal.II/include/deal.II/lac/relaxation_block.templates.h
deal.II/include/deal.II/lac/trilinos_block_vector.h
deal.II/include/deal.II/lac/trilinos_parallel_block_vector.h
deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
deal.II/include/deal.II/lac/trilinos_vector.h
deal.II/include/deal.II/lac/trilinos_vector_base.h
deal.II/include/deal.II/matrix_free/matrix_free.h
deal.II/include/deal.II/meshworker/dof_info.h
deal.II/include/deal.II/meshworker/integration_info.h
deal.II/include/deal.II/meshworker/loop.h
deal.II/include/deal.II/meshworker/simple.h
deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
deal.II/include/deal.II/multigrid/mg_tools.h
deal.II/include/deal.II/numerics/vector_tools.templates.h
deal.II/source/lac/block_sparsity_pattern.cc
deal.II/source/lac/trilinos_block_vector.cc
deal.II/source/lac/trilinos_sparse_matrix.cc
deal.II/source/lac/trilinos_sparsity_pattern.cc
deal.II/source/lac/trilinos_vector.cc
deal.II/source/lac/trilinos_vector_base.cc
deal.II/source/multigrid/mg_tools.cc
deal.II/source/multigrid/mg_tools.inst.in
deal.II/source/multigrid/mg_transfer_block.cc

index c2ba25e5539f2af06362026bffa007a231395aaa..b4ff8c8bf6e05e0267b62f7fd68308de8b6b1b85 100644 (file)
@@ -33,7 +33,12 @@ IF(FEATURE_BOOST_BUNDLED_CONFIGURED)
     )
 
   ADD_SUBDIRECTORY(${BOOST_FOLDER}/libs/serialization/src)
-  ADD_SUBDIRECTORY(${BOOST_FOLDER}/libs/iostreams/src)
+
+  IF(DEAL_II_WITH_ZLIB)
+    ADD_SUBDIRECTORY(${BOOST_FOLDER}/libs/iostreams/src)
+  ELSE()
+    MESSAGE(STATUS "BOOST::Iostreams will not be available because zlib or its header files could not be found")
+  ENDIF()
 
   IF(DEAL_II_WITH_THREADS AND NOT DEAL_II_WITH_CXX11)
     ADD_SUBDIRECTORY(${BOOST_FOLDER}/libs/thread/src)
index 07e336c8cca8e0de559df0f404e2855896147125..ef8f84838ce349ad5cf61a70e00073d49fbe9514 100644 (file)
 
 SET(src_boost_iostreams
     file_descriptor.cpp
-    gzip.cpp 
+    gzip.cpp
     mapped_file.cpp
-  )
-
-IF(DEAL_II_WITH_ZLIB)
-  SET(src_boost_iostreams
-      ${src_boost_iostreams}
-      zlib.cpp
-  )
-ELSE()
-  MESSAGE(STATUS "BOOST::Iostreams will not support gzipped streams because libz or its header files could not be found")
-ENDIF()
+    zlib.cpp
+    )
 
 IF(DEALII_WITH_BZIP2)
-  SET(src_boost_iostreams
-      ${src_boost_iostreams}
-      bzip2.cpp 
-  )
+  LIST(APPEND src_boost_iostreams bzip2.cpp)
 ELSE()
   MESSAGE(STATUS "BOOST::Iostreams will not support bz2'ed streams because libbz2 or its header files could not be found")
 ENDIF()
 
-
 DEAL_II_ADD_LIBRARY(obj_boost_iostreams OBJECT ${src_boost_iostreams})
-
index 903c19235df796259a0cde257aa964ddfa470d08..378851f71cfb5c35635711555746d343fba478da 100644 (file)
@@ -51,8 +51,11 @@ SET(CMAKE_PREFIX_PATH ${_cmake_prefix_path_backup})
 
 MARK_AS_ADVANCED(
   atlas_LIBRARY
+  atlcblas_LIBRARY
+  atllapack_LIBRARY
   blas_LIBRARY
   eigen_blas_LIBRARY
+  f77blas_LIBRARY
   gslcblas_LIBRARY
   lapack_LIBRARY
   m_LIBRARY
index 3cad88e7d21be58ea0e3d7f25d57eb6349e3ae8a..c4fc91f1dd986fc0254585b78cf23a998c638747 100644 (file)
 #
 #     DEAL_II_WITH_64BIT_INDICES
 #
+# Miscellaneous options:
+#     DEAL_II_DOXYGEN_USE_MATHJAX
+#
+#
 # *)  May also be set via environment variable (CXXFLAGS, LDFLAGS)
 #     (a nonempty cached variable has precedence and will not be
 #     overwritten by environment)
@@ -331,6 +335,12 @@ OPTION(DEAL_II_WITH_64BIT_INDICES
 #                                                                      #
 ########################################################################
 
+OPTION(DEAL_II_DOXYGEN_USE_MATHJAX
+  "If set to ON, doxygen documentation is generated using mathjax"
+  OFF
+  )
+MARK_AS_ADVANCED(DEAL_II_DOXYGEN_USE_MATHJAX)
+
 #
 # We do not support installation into the binary directory any more ("too
 # much pain, not enough profit"):
index cdeeeebb782ef3d6a51ccb8ba962752da25ec338..e813105703504d08c5cc13f522f8efc4c9cd8a29 100644 (file)
@@ -50,6 +50,11 @@ ADD_SUBDIRECTORY(tutorial) # uses ${DEAL_II_STEPS}
 #
 # Prepare auxiliary files for doxygen:
 #
+IF(DEAL_II_DOXYGEN_USE_MATHJAX)
+  SET(_use_mathjax YES)
+ELSE()
+  SET(_use_mathjax NO)
+ENDIF()
 CONFIGURE_FILE(
   ${CMAKE_CURRENT_SOURCE_DIR}/options.dox.in
   ${CMAKE_CURRENT_BINARY_DIR}/options.dox
index 2edcd2d5a8b85dfc5698642cd49ba917353164cf..b2ad033baa1e5d5d8601f1e131c6a56af2d6a6f1 100644 (file)
  * <b>v</b>=(0,0,<i>v</i><sub>3</sub>)<sup>T</sup>. writing the outcomes below
  * each other, we obtain the system
  * @f[
- * \arraycolsep1pt
  * \begin{matrix}
  * (\nabla u_1,\nabla v_1) &&& = (f_1, v_1)
  * \\
  * exhibit this system structure. Let us close by writing the full
  * system of the elastic equation with symmetric gradient <b>D</b>:
  * @f[
- * \arraycolsep1pt
  * \begin{matrix}
  * (\nabla u_1,\nabla v_1) + (\partial_1 u_1,\partial_1 v_1)
  * & (\partial_1 u_2,\partial_2 v_1)
index 851d8954c8283e5103e297bb41d78677418e22e0..88ca8938d2f05e240318bc96bb841910cf5f43ce 100644 (file)
@@ -98,7 +98,7 @@ HTML_OUTPUT            = deal.II
 HTML_FILE_EXTENSION    = .html
 HTML_HEADER            = header.html
 HTML_FOOTER            = footer.html
-USE_MATHJAX            = NO
+USE_MATHJAX            = @_use_mathjax@
 MATHJAX_RELPATH        = http://cdn.mathjax.org/mathjax/latest
 MATHJAX_EXTENSIONS     = TeX/AMSmath TeX/AMSsymbols
 HTML_COLORSTYLE_HUE    = 220
index e398c01f87cf74c1fc2d2d8e9491a2a389a863e7..8b705acd64638047a5739a5365ef347914720c57 100644 (file)
@@ -1,4 +1,4 @@
-// ---------------------------------------------------------------------
+ // ---------------------------------------------------------------------
 // $Id$
 //
 // Copyright (C) 2013 by the deal.II authors
@@ -106,7 +106,9 @@ inconvenience this causes.
   written into as an optional argument. This allows for pre-allocating all
   possible entries right away, which makes writing into the matrix from
   several threads possible (otherwise, only one processor at a time can write
-  off-processor data).
+  off-processor data). Similarly, TrilinosWrappers::MPI::Vector objects can
+  be initialized with hints to ghost elements for a writable vector that can
+  be added into from multiple threads.
   <br>
   (Martin Kronbichler, 2013/12/23)
   </li>
index c6a93e4c2ddda2d7d41411f1e3ab1f04d1b8b1fc..f4d43108b497656da05bad32c3abe6fe85460cfa 100644 (file)
 }
 
 
+@article{dealII81,
+  title = {The {\tt deal.{I}{I}} Library, Version 8.1},
+  author = {W. Bangerth and T. Heister and L. Heltai and G. Kanschat
+   and M. Kronbichler and M. Maier and B. Turcksin and T. D. Young},
+  journal = {arXiv preprint \url{http://arxiv.org/abs/1312.2266v4}},
+  year = {2013},
+}
+
 @article{dealII80,
   title = {The {\tt deal.{I}{I}} Library, Version 8.0},
   author = {W. Bangerth and T. Heister and L. Heltai and G. Kanschat
    and M. Kronbichler and M. Maier and B. Turcksin and T. D. Young},
-  journal = {arXiv preprint \url{http://arxiv.org/abs/1312.2266}},
+  journal = {arXiv preprint \url{http://arxiv.org/abs/1312.2266v3}},
   year = {2013},
 }
 
index 1c495eb61cdc3f3122b750d6892403db28dc44ee..cbcfc738746dc44e280ff4968a3846f8e3277f15 100644 (file)
     </p>
 
     <ol>
+      <li>
+       Wolfgang Bangerth, Timo Heister, Luca Heltai, Guido Kanschat, Martin Kronbichler, Matthias Maier, and Toby D. Young <br/>
+        <strong>The <acronym>deal.II</acronym> Library, Version 8.1
+        </strong>
+        <br>
+       <a href="http://arxiv.org/abs/1312.2266v4" target="_top">arxiv:1312.2266v4</a>, 2013.
+       <br>
+<pre>
+@article{dealII81,
+  title = {The {\tt deal.{I}{I}} Library, Version 8.1},
+  author = {W. Bangerth and T. Heister and L. Heltai and G. Kanschat
+   and M. Kronbichler and M. Maier and B. Turcksin and T. D. Young},
+  journal = {arXiv preprint \url{http://arxiv.org/abs/1312.2266v4}},
+  year = {2013},
+}
+</pre>
+      </li>
       <li> <a href="http://www.math.tamu.edu/~bangerth"
               target="_top">W. Bangerth</a>,
        <a href="http://ganymed.iwr.uni-heidelberg.de/~hartmann"
 }
 </pre>
       </li>
+    </ol>
+
+    <p>
+      For your convenience, we provide a <a href="dealii.bib">BibTeX file</a> with entries
+      for the publications in this section.
+    </p>
+
+
+    <h4>Older Releases</h4>
+    
+    <ol>
       <li>
        Wolfgang Bangerth, Timo Heister, Luca Heltai, Guido Kanschat, Martin Kronbichler, Matthias Maier, and Toby D. Young <br/>
         <strong>The <acronym>deal.II</acronym> Library, Version 8.0
         </strong>
         <br>
-       <a href="http://arxiv.org/abs/1312.2266" target="_top">arxiv:1312.2266</a>, 2013.
+       <a href="http://arxiv.org/abs/1312.2266v3" target="_top">arxiv:1312.2266v3</a>, 2013.
        <br>
 <pre>
 @article{dealII80,
   title = {The {\tt deal.{I}{I}} Library, Version 8.0},
   author = {W. Bangerth and T. Heister and L. Heltai and G. Kanschat
    and M. Kronbichler and M. Maier and B. Turcksin and T. D. Young},
-  journal = {arXiv preprint \url{http://arxiv.org/abs/1312.2266}},
+  journal = {arXiv preprint \url{http://arxiv.org/abs/1312.2266v3}},
   year = {2013},
 }
 </pre>
       </li>
     </ol>
 
-    <p>
-      For your convenience, we provide a <a href="dealii.bib">BibTeX file</a> with entries
-      for the publications in this section.
-    </p>
-
-
 
     <!-- check on Freddi papers, conference at http://www.eurosim2007.org -->
     <!-- watch for the papers here -->
index 26c26d3acf88092228e297ce4113b382ffbc6175..1bae0ad4a543ac705dcf73ee110071268263d1bf 100644 (file)
@@ -306,7 +306,7 @@ namespace Step12
      &AdvectionProblem<dim>::integrate_cell_term,
      &AdvectionProblem<dim>::integrate_boundary_term,
      &AdvectionProblem<dim>::integrate_face_term,
-     assembler, true);
+     assembler);
   }
 
 
index e22741f97c971659a2df15d6a874524b901a8ea7..8082826f243a580dfbf5b9a2ae6255c13d19225b 100644 (file)
@@ -1060,6 +1060,22 @@ namespace TrilinosWrappers
     BlockSparsityPattern (const std::vector<IndexSet> &parallel_partitioning,
                           const MPI_Comm &communicator = MPI_COMM_WORLD);
 
+    /**
+     * Initialize the pattern with two arrays of index sets that specify rows
+     * and columns of the matrix, where the size() of the IndexSets specifies
+     * the size of the blocks and the values in each IndexSet denotes the rows
+     * that are going to be saved in each block. The additional index set
+     * writable_rows is used to set all rows that we allow to write
+     * locally. This constructor is used to create matrices that allow several
+     * threads to write simultaneously into the matrix (to different rows, of
+     * course), see the method TrilinosWrappers::SparsityPattern::reinit
+     * method with three index set arguments for more details.
+     */
+    BlockSparsityPattern (const std::vector<IndexSet> &row_parallel_partitioning,
+                          const std::vector<IndexSet> &column_parallel_partitioning,
+                          const std::vector<IndexSet> &writeable_rows,
+                          const MPI_Comm              &communicator = MPI_COMM_WORLD);
+
     /**
      * Resize the matrix to a tensor product of matrices with dimensions
      * defined by the arguments.
index 6e6bb9bb159517be70a041509f55a7d36633718a..994a5c502c7c43d414c3575e4884425a6f751a31 100644 (file)
@@ -86,19 +86,6 @@ public:
                     const bool invert_diagonal = true,
                     const bool same_diagonal = false);
 
-    /**
-     * @deprecated Since the BlockList is now a data member of this
-     * class, it is not recommended anymore to generate it
-     * independently.
-     *
-     * Constructor.
-     */
-    AdditionalData (const BlockList &block_list,
-                    const double relaxation = 1.,
-                    const bool invert_diagonal = true,
-                    const bool same_diagonal = false) DEAL_II_DEPRECATED;
-
-
     /**
      * The mapping from indices to blocks.
      */
index 621ed762806e2f2ec468d0aa92a5b91eaac706ec..8028f74acf6de7b83769cd90ad730d19ad4d3c8e 100644 (file)
@@ -51,24 +51,6 @@ RelaxationBlock<MATRIX,inverse_type>::AdditionalData::memory_consumption() const
 }
 
 
-template <class MATRIX, typename inverse_type>
-inline
-RelaxationBlock<MATRIX,inverse_type>::AdditionalData::AdditionalData (
-  const BlockList &bl,
-  const double relaxation,
-  const bool invert_diagonal,
-  const bool same_diagonal)
-  :
-  relaxation(relaxation),
-  invert_diagonal(invert_diagonal),
-  same_diagonal(same_diagonal),
-  inversion(PreconditionBlockBase<inverse_type>::gauss_jordan),
-  threshold(0.)
-{
-  bl.create_sparsity_pattern(block_list, 0);
-}
-
-
 template <class MATRIX, typename inverse_type>
 inline
 void
index 9006b134c22d983fb62b9fd8dc8c44de4dc0c840..276f778902637c1292d34b091a64ad272fffcbca 100644 (file)
@@ -71,20 +71,17 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Typedef the base class for simpler
-     * access to its own typedefs.
+     * Typedef the base class for simpler access to its own typedefs.
      */
     typedef BlockVectorBase<Vector> BaseClass;
 
     /**
-     * Typedef the type of the underlying
-     * vector.
+     * Typedef the type of the underlying vector.
      */
     typedef BaseClass::BlockType  BlockType;
 
     /**
-     * Import the typedefs from the base
-     * class.
+     * Import the typedefs from the base class.
      */
     typedef BaseClass::value_type      value_type;
     typedef BaseClass::pointer         pointer;
@@ -96,90 +93,59 @@ namespace TrilinosWrappers
     typedef BaseClass::const_iterator  const_iterator;
 
     /**
-     * Default constructor. Generate an
-     * empty vector without any blocks.
+     * Default constructor. Generate an empty vector without any blocks.
      */
     BlockVector ();
 
     /**
-     * Constructor. Generate a block
-     * vector with as many blocks as
-     * there are entries in Input_Maps.
-     * For this non-distributed vector,
-     * the %parallel partitioning is not
-     * used, just the global size of the
-     * partitioner.
+     * Constructor. Generate a block vector with as many blocks as there are
+     * entries in Input_Maps.  For this non-distributed vector, the %parallel
+     * partitioning is not used, just the global size of the partitioner.
      */
     explicit BlockVector (const std::vector<Epetra_Map> &partitioner);
 
     /**
-     * Constructor. Generate a block
-     * vector with as many blocks as
-     * there are entries in Input_Maps.
-     * For this non-distributed vector,
-     * the %parallel partitioning is not
-     * used, just the global size of the
-     * partitioner.
+     * Constructor. Generate a block vector with as many blocks as there are
+     * entries in Input_Maps.  For this non-distributed vector, the %parallel
+     * partitioning is not used, just the global size of the partitioner.
      */
     explicit BlockVector (const std::vector<IndexSet> &partitioner,
                           const MPI_Comm              &communicator = MPI_COMM_WORLD);
 
     /**
-     * Copy-Constructor. Set all the
-     * properties of the non-%parallel
-     * vector to those of the given
-     * %parallel vector and import the
-     * elements.
+     * Copy-Constructor. Set all the properties of the non-%parallel vector to
+     * those of the given %parallel vector and import the elements.
      */
     BlockVector (const MPI::BlockVector &V);
 
     /**
-     * Copy-Constructor. Set all the
-     * properties of the vector to those
-     * of the given input vector and copy
-     * the elements.
+     * Copy-Constructor. Set all the properties of the vector to those of the
+     * given input vector and copy the elements.
      */
     BlockVector (const BlockVector  &V);
 
     /**
-     * Creates a block vector
-     * consisting of
-     * <tt>num_blocks</tt>
-     * components, but there is no
-     * content in the individual
-     * components and the user has to
-     * fill appropriate data using a
-     * reinit of the blocks.
+     * Creates a block vector consisting of <tt>num_blocks</tt> components,
+     * but there is no content in the individual components and the user has
+     * to fill appropriate data using a reinit of the blocks.
      */
     explicit BlockVector (const size_type num_blocks);
 
     /**
-     * Constructor. Set the number of
-     * blocks to <tt>n.size()</tt> and
-     * initialize each block with
-     * <tt>n[i]</tt> zero elements.
+     * Constructor. Set the number of blocks to <tt>n.size()</tt> and
+     * initialize each block with <tt>n[i]</tt> zero elements.
      *
      * References BlockVector.reinit().
      */
     explicit BlockVector (const std::vector<size_type> &N);
 
     /**
-     * Constructor. Set the number of
-     * blocks to
-     * <tt>n.size()</tt>. Initialize the
-     * vector with the elements
-     * pointed to by the range of
-     * iterators given as second and
-     * third argument. Apart from the
-     * first argument, this
-     * constructor is in complete
-     * analogy to the respective
-     * constructor of the
-     * <tt>std::vector</tt> class, but the
-     * first argument is needed in
-     * order to know how to subdivide
-     * the block vector into
-     * different blocks.
+     * Constructor. Set the number of blocks to <tt>n.size()</tt>. Initialize
+     * the vector with the elements pointed to by the range of iterators given
+     * as second and third argument. Apart from the first argument, this
+     * constructor is in complete analogy to the respective constructor of the
+     * <tt>std::vector</tt> class, but the first argument is needed in order
+     * to know how to subdivide the block vector into different blocks.
      */
     template <typename InputIterator>
     BlockVector (const std::vector<size_type> &n,
@@ -196,8 +162,7 @@ namespace TrilinosWrappers
     *
     * @deprecated
     *
-    * See @ref GlossCompress "Compressing
-    * distributed objects" for more
+    * See @ref GlossCompress "Compressing distributed objects" for more
     * information.
     */
     void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED;
@@ -208,91 +173,60 @@ namespace TrilinosWrappers
     using BlockVectorBase<Vector>::compress;
 
     /**
-     * Copy operator: fill all
-     * components of the vector that
-     * are locally stored with the
-     * given scalar value.
+     * Copy operator: fill all components of the vector that are locally
+     * stored with the given scalar value.
      */
     BlockVector &
     operator = (const value_type s);
 
     /**
-     * Copy operator for a
-     * distributed Trilinos vector to
-     * a localized one.
+     * Copy operator for a distributed Trilinos vector to a localized one.
      */
     BlockVector &
     operator = (const MPI::BlockVector &V);
 
     /**
-     * Copy operator for arguments of
-     * the same type.
+     * Copy operator for arguments of the same type.
      */
     BlockVector &
     operator = (const BlockVector &V);
 
     /**
-     * Another copy function. This
-     * one takes a deal.II block
-     * vector and copies it into a
-     * TrilinosWrappers block
-     * vector. Note that the number
-     * of blocks has to be the same
-     * in the vector as in the input
-     * vector. Use the reinit()
-     * command for resizing the
-     * BlockVector or for changing
-     * the internal structure of the
-     * block components.
+     * Another copy function. This one takes a deal.II block vector and copies
+     * it into a TrilinosWrappers block vector. Note that the number of blocks
+     * has to be the same in the vector as in the input vector. Use the
+     * reinit() command for resizing the BlockVector or for changing the
+     * internal structure of the block components.
      *
-     * Since Trilinos only works on
-     * doubles, this function is
-     * limited to accept only one
-     * possible number type in the
-     * deal.II vector.
+     * Since Trilinos only works on doubles, this function is limited to
+     * accept only one possible number type in the deal.II vector.
      */
     template <typename Number>
     BlockVector &
     operator = (const ::dealii::BlockVector<Number> &V);
 
     /**
-     * Reinitialize the BlockVector to
-     * contain as many blocks as there
-     * are Epetra_Maps given in the
-     * input argument, according to the
-     * global size of the individual
-     * components described in the
-     * maps. Note that the resulting
-     * vector will be stored completely
-     * on each process. The Epetra_Map
-     * is useful when data exchange
-     * with a distributed vector based
-     * on the same Epetra_map is
-     * intended. In that case, the same
-     * communicator is used for data
-     * exchange.
+     * Reinitialize the BlockVector to contain as many blocks as there are
+     * Epetra_Maps given in the input argument, according to the global size
+     * of the individual components described in the maps. Note that the
+     * resulting vector will be stored completely on each process. The
+     * Epetra_Map is useful when data exchange with a distributed vector based
+     * on the same Epetra_map is intended. In that case, the same communicator
+     * is used for data exchange.
      *
-     * If <tt>fast==false</tt>, the vector
-     * is filled with zeros.
+     * If <tt>fast==false</tt>, the vector is filled with zeros.
      */
     void reinit (const std::vector<Epetra_Map> &partitioning,
                  const bool                     fast = false);
 
     /**
-     * Reinitialize the BlockVector to
-     * contain as many blocks as there
-     * are index sets given in the
-     * input argument, according to the
-     * global size of the individual
-     * components described in the
-     * index set, and using a given MPI
-     * communicator. The MPI
-     * communicator is useful when data
-     * exchange with a distributed
-     * vector based on the same
-     * initialization is intended. In
-     * that case, the same communicator
-     * is used for data exchange.
+     * Reinitialize the BlockVector to contain as many blocks as there are
+     * index sets given in the input argument, according to the global size of
+     * the individual components described in the index set, and using a given
+     * MPI communicator. The MPI communicator is useful when data exchange
+     * with a distributed vector based on the same initialization is
+     * intended. In that case, the same communicator is used for data
+     * exchange.
      *
      * If <tt>fast==false</tt>, the vector
      * is filled with zeros.
@@ -302,13 +236,9 @@ namespace TrilinosWrappers
                  const bool                   fast = false);
 
     /**
-     * Reinitialize the BlockVector to
-     * contain as many blocks as there
-     * are elements in the first
-     * argument, and with the respective
-     * sizes. Since no distribution map
-     * is given, all vectors are local
-     * vectors.
+     * Reinitialize the BlockVector to contain as many blocks as there are
+     * elements in the first argument, and with the respective sizes. Since no
+     * distribution map is given, all vectors are local vectors.
      *
      * If <tt>fast==false</tt>, the vector
      * is filled with zeros.
@@ -324,76 +254,44 @@ namespace TrilinosWrappers
     void reinit (const MPI::BlockVector &V);
 
     /**
-     * Change the dimension to that
-     * of the vector <tt>V</tt>. The same
-     * applies as for the other
-     * reinit() function.
+     * Change the dimension to that of the vector <tt>V</tt>. The same applies
+     * as for the other reinit() function.
      *
-     * The elements of <tt>V</tt> are not
-     * copied, i.e.  this function is
-     * the same as calling <tt>reinit
-     * (V.size(), fast)</tt>.
+     * The elements of <tt>V</tt> are not copied, i.e.  this function is the
+     * same as calling <tt>reinit (V.size(), fast)</tt>.
      *
-     * Note that you must call this
-     * (or the other reinit()
-     * functions) function, rather
-     * than calling the reinit()
-     * functions of an individual
-     * block, to allow the block
-     * vector to update its caches of
-     * vector sizes. If you call
-     * reinit() on one of the
-     * blocks, then subsequent
-     * actions on this object may
-     * yield unpredictable results
-     * since they may be routed to
-     * the wrong block.
+     * Note that you must call this (or the other reinit() functions)
+     * function, rather than calling the reinit() functions of an individual
+     * block, to allow the block vector to update its caches of vector
+     * sizes. If you call reinit() on one of the blocks, then subsequent
+     * actions on this object may yield unpredictable results since they may
+     * be routed to the wrong block.
      */
     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.
+     * 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 size_type num_blocks);
 
     /**
-     * Swap the contents of this
-     * vector and the other vector
-     * <tt>v</tt>. One could do this
-     * operation with a temporary
-     * variable and copying over the
-     * data elements, but this
-     * function is significantly more
-     * efficient since it only swaps
-     * the pointers to the data of
-     * the two vectors and therefore
-     * does not need to allocate
-     * temporary storage and move
-     * data around.
+     * Swap the contents of this vector and the other vector <tt>v</tt>. One
+     * could do this operation with a temporary variable and copying over the
+     * data elements, but this function is significantly more efficient since
+     * it only swaps the pointers to the data of the two vectors and therefore
+     * does not need to allocate temporary storage and move data around.
      *
-     * Limitation: right now this
-     * function only works if both
-     * vectors have the same number
-     * of blocks. If needed, the
-     * numbers of blocks should be
+     * Limitation: right now this function only works if both vectors have the
+     * same number of blocks. If needed, the numbers of blocks should be
      * exchanged, too.
      *
-     * This function is analog to the
-     * the swap() function of all C++
-     * standard containers. Also,
-     * there is a global function
-     * swap(u,v) that simply calls
-     * <tt>u.swap(v)</tt>, again in analogy
-     * to standard functions.
+     * This function is analog to the the swap() function of all C++ standard
+     * containers. Also, there is a global function swap(u,v) that simply
+     * calls <tt>u.swap(v)</tt>, again in analogy to standard functions.
      */
     void swap (BlockVector &v);
 
index 7ec7a0e217891eac190dc548609ee3366ccd4f63..4e990c1469a3a0b1470441ff6b05fef9c3b9a3ef 100644 (file)
@@ -73,20 +73,17 @@ namespace TrilinosWrappers
     {
     public:
       /**
-       * Typedef the base class for simpler
-       * access to its own typedefs.
+       * Typedef the base class for simpler access to its own typedefs.
        */
       typedef BlockVectorBase<Vector> BaseClass;
 
       /**
-       * Typedef the type of the underlying
-       * vector.
+       * Typedef the type of the underlying vector.
        */
       typedef BaseClass::BlockType  BlockType;
 
       /**
-       * Import the typedefs from the base
-       * class.
+       * Import the typedefs from the base class.
        */
       typedef BaseClass::value_type      value_type;
       typedef BaseClass::pointer         pointer;
@@ -98,62 +95,46 @@ namespace TrilinosWrappers
       typedef BaseClass::const_iterator  const_iterator;
 
       /**
-       * Default constructor. Generate an
-       * empty vector without any blocks.
+       * Default constructor. Generate an empty vector without any blocks.
        */
       BlockVector ();
 
       /**
-       * Constructor. Generate a block
-       * vector with as many blocks as
-       * there are entries in @p
-       * partitioning. Each Epetra_Map
-       * contains the layout of the
-       * distribution of data among the MPI
-       * processes.
+       * Constructor. Generate a block vector with as many blocks as there are
+       * entries in @p partitioning. Each Epetra_Map contains the layout of
+       * the distribution of data among the MPI processes.
        */
       explicit BlockVector (const std::vector<Epetra_Map> &parallel_partitioning);
 
       /**
-       * Constructor. Generate a block
-       * vector with as many blocks as
-       * there are entries in
-       * @p partitioning.  Each IndexSet
-       * together with the MPI communicator
-       * contains the layout of the
-       * distribution of data among the MPI
-       * processes.
+       * Constructor. Generate a block vector with as many blocks as there are
+       * entries in @p partitioning.  Each IndexSet together with the MPI
+       * communicator contains the layout of the distribution of data among
+       * the MPI processes.
        */
       explicit BlockVector (const std::vector<IndexSet> &parallel_partitioning,
                             const MPI_Comm              &communicator = MPI_COMM_WORLD);
 
       /**
-       * Creates a BlockVector with ghost elements. @p ghost_values
-       * may contain any elements in @p parallel_partitioning, they will
-       * be ignored.
+       * Creates a BlockVector with ghost elements. See the respective
+       * reinit() method for more details. @p ghost_values may contain any
+       * elements in @p parallel_partitioning, they will be ignored.
        */
       BlockVector (const std::vector<IndexSet> &parallel_partitioning,
                    const std::vector<IndexSet> &ghost_values,
-                   const MPI_Comm              &communicator);
-
+                   const MPI_Comm              &communicator,
+                   const bool                   vector_writable = false);
 
       /**
-       * Copy-Constructor. Set all the
-       * properties of the parallel vector
-       * to those of the given argument and
-       * copy the elements.
+       * Copy-Constructor. Set all the properties of the parallel vector to
+       * those of the given argument and copy the elements.
        */
       BlockVector (const BlockVector  &V);
 
       /**
-       * Creates a block vector
-       * consisting of
-       * <tt>num_blocks</tt>
-       * components, but there is no
-       * content in the individual
-       * components and the user has to
-       * fill appropriate data using a
-       * reinit of the blocks.
+       * Creates a block vector consisting of <tt>num_blocks</tt> components,
+       * but there is no content in the individual components and the user has
+       * to fill appropriate data using a reinit of the blocks.
        */
       explicit BlockVector (const size_type num_blocks);
 
@@ -163,168 +144,124 @@ namespace TrilinosWrappers
       ~BlockVector ();
 
       /**
-       * Copy operator: fill all
-       * components of the vector that
-       * are locally stored with the
-       * given scalar value.
+       * Copy operator: fill all components of the vector that are locally
+       * stored with the given scalar value.
        */
       BlockVector &
       operator = (const value_type s);
 
       /**
-       * Copy operator for arguments of
-       * the same type.
+       * Copy operator for arguments of the same type.
        */
       BlockVector &
       operator = (const BlockVector &V);
 
       /**
-       * Copy operator for arguments of
-       * the localized Trilinos vector
-       * type.
+       * Copy operator for arguments of the localized Trilinos vector type.
        */
       BlockVector &
       operator = (const ::dealii::TrilinosWrappers::BlockVector &V);
 
       /**
-       * Another copy function. This
-       * one takes a deal.II block
-       * vector and copies it into a
-       * TrilinosWrappers block
-       * vector. Note that the number
-       * of blocks has to be the same
-       * in the vector as in the input
-       * vector. Use the reinit()
-       * command for resizing the
-       * BlockVector or for changing
-       * the internal structure of the
-       * block components.
+       * Another copy function. This one takes a deal.II block vector and
+       * copies it into a TrilinosWrappers block vector. Note that the number
+       * of blocks has to be the same in the vector as in the input
+       * vector. Use the reinit() command for resizing the BlockVector or for
+       * changing the internal structure of the block components.
        *
-       * Since Trilinos only works on
-       * doubles, this function is
-       * limited to accept only one
-       * possible number type in the
-       * deal.II vector.
+       * Since Trilinos only works on doubles, this function is limited to
+       * accept only one possible number type in the deal.II vector.
        */
       template <typename Number>
       BlockVector &
       operator = (const ::dealii::BlockVector<Number> &V);
 
       /**
-       * Reinitialize the BlockVector to
-       * contain as many blocks as there
-       * are Epetra_Maps given in the input
-       * argument, according to the
-       * parallel distribution of the
-       * individual components described
-       * in the maps.
+       * Reinitialize the BlockVector to contain as many blocks as there are
+       * Epetra_Maps given in the input argument, according to the parallel
+       * distribution of the individual components described in the maps.
        *
-       * If <tt>fast==false</tt>, the vector
-       * is filled with zeros.
+       * If <tt>fast==false</tt>, the vector is filled with zeros.
        */
       void reinit (const std::vector<Epetra_Map> &parallel_partitioning,
                    const bool                     fast = false);
 
       /**
-       * Reinitialize the BlockVector to
-       * contain as many blocks as there
-       * are index sets given in the input
-       * argument, according to the
-       * parallel distribution of the
-       * individual components described
-       * in the maps.
+       * Reinitialize the BlockVector to contain as many blocks as there are
+       * index sets given in the input argument, according to the parallel
+       * distribution of the individual components described in the maps.
        *
-       * If <tt>fast==false</tt>, the vector
-       * is filled with zeros.
+       * If <tt>fast==false</tt>, the vector is filled with zeros.
        */
       void reinit (const std::vector<IndexSet> &parallel_partitioning,
                    const MPI_Comm              &communicator = MPI_COMM_WORLD,
                    const bool                   fast = false);
+
       /**
-       * like above, but with a second set of indices for
-       * ghost entries.
+       * Reinit functionality. This function destroys the old vector content
+       * and generates a new one based on the input partitioning. In addition
+       * to just specifying one index set as in all the other methods above,
+       * this method allows to supply an additional set of ghost
+       * entries. There are two different versions of a vector that can be
+       * created. If the flag @p vector_writable is set to @p false, the
+       * vector only allows read access to the joint set of @p
+       * parallel_partitioning and @p ghost_entries. The effect of the reinit
+       * method is then equivalent to calling the other reinit method with an
+       * index set containing both the locally owned entries and the ghost
+       * entries.
+       *
+       * If the flag @p vector_writable is set to true, this creates an
+       * alternative storage scheme for ghost elements that allows multiple
+       * threads to write into the vector (for the other reinit methods, only
+       * one thread is allowed to write into the ghost entries at a time).
        */
       void reinit (const std::vector<IndexSet> &partitioning,
                    const std::vector<IndexSet> &ghost_values,
-                   const MPI_Comm              &communicator = MPI_COMM_WORLD);
+                   const MPI_Comm              &communicator = MPI_COMM_WORLD,
+                   const bool      vector_writable = false);
 
 
       /**
-       * Change the dimension to that
-       * of the vector <tt>V</tt>. The same
-       * applies as for the other
-       * reinit() function.
+       * Change the dimension to that of the vector <tt>V</tt>. The same
+       * applies as for the other reinit() function.
        *
-       * The elements of <tt>V</tt> are not
-       * copied, i.e.  this function is
-       * the same as calling <tt>reinit
-       * (V.size(), fast)</tt>.
+       * The elements of <tt>V</tt> are not copied, i.e.  this function is the
+       * same as calling <tt>reinit (V.size(), fast)</tt>.
        *
-       * Note that you must call this
-       * (or the other reinit()
-       * functions) function, rather
-       * than calling the reinit()
-       * functions of an individual
-       * block, to allow the block
-       * vector to update its caches of
-       * vector sizes. If you call
-       * reinit() on one of the
-       * blocks, then subsequent
-       * actions on this object may
-       * yield unpredictable results
-       * since they may be routed to
-       * the wrong block.
+       * Note that you must call this (or the other reinit() functions)
+       * function, rather than calling the reinit() functions of an individual
+       * block, to allow the block vector to update its caches of vector
+       * sizes. If you call reinit() on one of the blocks, then subsequent
+       * actions on this object may yield unpredictable results since they may
+       * be routed to the wrong block.
        */
       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.
+       * 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 size_type num_blocks);
 
       /**
-       * This reinit function is meant to
-       * be used for parallel
-       * calculations where some
-       * non-local data has to be
-       * used. The typical situation
-       * where one needs this function is
-       * the call of the
-       * FEValues<dim>::get_function_values
-       * function (or of some
-       * derivatives) in parallel. Since
-       * it is usually faster to retrieve
-       * the data in advance, this
-       * function can be called before
-       * the assembly forks out to the
-       * different processors. What this
-       * function does is the following:
-       * It takes the information in the
-       * columns of the given matrix and
-       * looks which data couples between
-       * the different processors. That
-       * data is then queried from the
-       * input vector. Note that you
-       * should not write to the
-       * resulting vector any more, since
-       * the some data can be stored
-       * several times on different
-       * processors, leading to
-       * unpredictable results. In
-       * particular, such a vector cannot
-       * be used for matrix-vector
-       * products as for example done
-       * during the solution of linear
-       * systems.
+       * This reinit function is meant to be used for parallel calculations
+       * where some non-local data has to be used. The typical situation where
+       * one needs this function is the call of the
+       * FEValues<dim>::get_function_values function (or of some derivatives)
+       * in parallel. Since it is usually faster to retrieve the data in
+       * advance, this function can be called before the assembly forks out to
+       * the different processors. What this function does is the following:
+       * It takes the information in the columns of the given matrix and looks
+       * which data couples between the different processors. That data is
+       * then queried from the input vector. Note that you should not write to
+       * the resulting vector any more, since the some data can be stored
+       * several times on different processors, leading to unpredictable
+       * results. In particular, such a vector cannot be used for
+       * matrix-vector products as for example done during the solution of
+       * linear systems.
        */
       void import_nonlocal_data_for_fe (const TrilinosWrappers::BlockSparseMatrix &m,
                                         const BlockVector                         &v);
@@ -335,8 +272,7 @@ namespace TrilinosWrappers
        *
        * @deprecated
        *
-       * See @ref GlossCompress "Compressing
-       * distributed objects" for more
+       * See @ref GlossCompress "Compressing distributed objects" for more
        * information.
        */
       void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED;
@@ -348,17 +284,10 @@ namespace TrilinosWrappers
 
 
       /**
-       * Returns the state of the
-       * vector, i.e., whether
-       * compress() needs to be
-       * called after an operation
-       * requiring data
-       * exchange. Does only return
-       * non-true values when used in
-       * <tt>debug</tt> mode, since
-       * it is quite expensive to
-       * keep track of all operations
-       * that lead to the need for
+       * Returns the state of the vector, i.e., whether compress() needs to be
+       * called after an operation requiring data exchange. Does only return
+       * non-true values when used in <tt>debug</tt> mode, since it is quite
+       * expensive to keep track of all operations that lead to the need for
        * compress().
        */
       bool is_compressed () const;
@@ -369,34 +298,21 @@ namespace TrilinosWrappers
       bool has_ghost_elements() const;
 
       /**
-       * Swap the contents of this
-       * vector and the other vector
-       * <tt>v</tt>. One could do this
-       * operation with a temporary
-       * variable and copying over the
-       * data elements, but this
-       * function is significantly more
-       * efficient since it only swaps
-       * the pointers to the data of
-       * the two vectors and therefore
-       * does not need to allocate
-       * temporary storage and move
-       * data around.
+       * Swap the contents of this vector and the other vector <tt>v</tt>. One
+       * could do this operation with a temporary variable and copying over
+       * the data elements, but this function is significantly more efficient
+       * since it only swaps the pointers to the data of the two vectors and
+       * therefore does not need to allocate temporary storage and move data
+       * around.
        *
-       * Limitation: right now this
-       * function only works if both
-       * vectors have the same number
-       * of blocks. If needed, the
-       * numbers of blocks should be
+       * Limitation: right now this function only works if both vectors have
+       * the same number of blocks. If needed, the numbers of blocks should be
        * exchanged, too.
        *
-       * This function is analog to the
-       * the swap() function of all C++
-       * standard containers. Also,
-       * there is a global function
-       * swap(u,v) that simply calls
-       * <tt>u.swap(v)</tt>, again in analogy
-       * to standard functions.
+       * This function is analog to the the swap() function of all C++
+       * standard containers. Also, there is a global function swap(u,v) that
+       * simply calls <tt>u.swap(v)</tt>, again in analogy to standard
+       * functions.
        */
       void swap (BlockVector &v);
 
@@ -446,15 +362,19 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     BlockVector::BlockVector (const std::vector<IndexSet> &parallel_partitioning,
                               const std::vector<IndexSet> &ghost_values,
-                              const MPI_Comm              &communicator)
+                              const MPI_Comm              &communicator,
+                              const bool                   vector_writable)
     {
-      reinit(parallel_partitioning, ghost_values, communicator);
+      reinit(parallel_partitioning, ghost_values, communicator,
+             vector_writable);
     }
 
 
+
     inline
     BlockVector::BlockVector (const size_type num_blocks)
     {
@@ -533,6 +453,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     bool
     BlockVector::has_ghost_elements() const
@@ -545,6 +466,8 @@ namespace TrilinosWrappers
       return ghosted;
     }
 
+
+
     inline
     void
     BlockVector::swap (BlockVector &v)
index 8babcea0032da0b64826169943ee87ce9833e834..2a21d5daa9d8c10508ca50b7dfe66cb4c83dab27 100644 (file)
@@ -452,23 +452,23 @@ namespace TrilinosWrappers
    *
    * When writing into Trilinos matrices from several threads in shared
    * memory, several things must be kept in mind as there is no built-in locks
-   * in this class to prevent data races. Therefore, simultaneous access to
-   * the same matrix row at the same time leads to data races in general and
-   * must be explicitly avoided by the user. However, it is possible to access
-   * <b>different</b> rows of the matrix from several threads simultaneously
-   * under the following two conditions:
+   * in this class to prevent data races. Simultaneous access to the same
+   * matrix row at the same time can lead to data races and must be explicitly
+   * avoided by the user. However, it is possible to access <b>different</b>
+   * rows of the matrix from several threads simultaneously under the
+   * following two conditions:
    * <ul>
    *   <li> The matrix uses only one MPI process.
    *   <li> The matrix has been initialized from a
    *   TrilinosWrappers::SparsityPattern object that in turn has been
    *   initialized with the reinit function specifying three index sets, one
    *   for the rows, one for the columns and for the larger set of @p
-   *   writeable_rows. Note that all other reinit methods and constructors of
-   *   TrilinosWrappers::SparsityPattern will result in a matrix that needs to
-   *   allocate off-processor entries on demand, which breaks
-   *   thread-safety. Of course, using the respective reinit method for the
-   *   block Trilinos sparsity pattern and block matrix also results in
-   *   thread-safety.
+   *   writeable_rows, and the operation is an addition. Note that all other
+   *   reinit methods and constructors of TrilinosWrappers::SparsityPattern
+   *   will result in a matrix that needs to allocate off-processor entries on
+   *   demand, which breaks thread-safety. Of course, using the respective
+   *   reinit method for the block Trilinos sparsity pattern and block matrix
+   *   also results in thread-safety.
    * </ul>
    *
    * @ingroup TrilinosWrappers
@@ -2499,25 +2499,17 @@ namespace TrilinosWrappers
             if (ierr > 0)
               ierr = 0;
           }
-        else if (nonlocal_matrix.get() != 0)
-          {
-            // this is the case when we have explicitly set the off-processor
-            // rows and want to create a separate matrix object for them (to
-            // retain thread-safety)
-            Assert (nonlocal_matrix->RowMap().LID(static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
-                    ExcMessage("Attempted to write into off-processor matrix row "
-                               "that has not be specified as being writable upon "
-                               "initialization"));
-            ierr = nonlocal_matrix->ReplaceGlobalValues(row, n_columns,
-                                                        col_value_ptr,
-                                                        col_index_ptr);
-          }
         else
           ierr = matrix->ReplaceGlobalValues (1,
                                               (TrilinosWrappers::types::int_type *)&row,
                                               n_columns, col_index_ptr,
                                               &col_value_ptr,
                                               Epetra_FECrsMatrix::ROW_MAJOR);
+        // use the FECrsMatrix facilities for set even in the case when we
+        // have explicitly set the off-processor rows because that only works
+        // properly when adding elements, not when setting them (since we want
+        // to only touch elements that have been set explicitly, and there is
+        // no way on the receiving processor to identify them otherwise)
       }
 
     Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
@@ -2762,8 +2754,7 @@ namespace TrilinosWrappers
 
 
 
-  // inline "simple" functions that are
-  // called frequently and do only involve
+  // inline "simple" functions that are called frequently and do only involve
   // a call to some Trilinos function.
   inline
   SparseMatrix::size_type
index 1c4762e18cd6a0e4b9eb063c96d073ef65bb7558..9e7896532b2066bfd985661cb04ae211693c1ee1 100644 (file)
@@ -186,6 +186,18 @@ namespace TrilinosWrappers
      * operations at the same time, for example by placing zero additions if
      * necessary.
      *
+     * <h3>Thread safety of Trilinos vectors</h3>
+     *
+     * When writing into Trilinos vectors from several threads in shared
+     * memory, several things must be kept in mind as there is no built-in
+     * locks in this class to prevent data races. Simultaneous access to the
+     * same vector entry at the same time results in data races and must be
+     * explicitly avoided by the user. However, it is possible to access
+     * <b>different</b> entries of the vector from several threads
+     * simultaneously when only one MPI process is present or the vector has
+     * been constructed with an additional index set for ghost entries in
+     * write mode.
+     *
      * @ingroup TrilinosWrappers
      * @ingroup Vectors
      * @author Martin Kronbichler, Wolfgang Bangerth, 2008, 2009
@@ -199,15 +211,14 @@ namespace TrilinosWrappers
       typedef dealii::types::global_dof_index size_type;
 
       /**
-       * A variable that indicates whether this vector
-       * supports distributed data storage. If true, then
-       * this vector also needs an appropriate compress()
-       * function that allows communicating recent set or
-       * add operations to individual elements to be communicated
-       * to other processors.
+       * A variable that indicates whether this vector supports distributed
+       * data storage. If true, then this vector also needs an appropriate
+       * compress() function that allows communicating recent set or add
+       * operations to individual elements to be communicated to other
+       * processors.
        *
-       * For the current class, the variable equals
-       * true, since it does support parallel data storage.
+       * For the current class, the variable equals true, since it does
+       * support parallel data storage.
        */
       static const bool supports_distributed_data = true;
 
@@ -216,20 +227,14 @@ namespace TrilinosWrappers
        */
       //@{
       /**
-       * Default constructor that
-       * generates an empty (zero size)
-       * vector. The function
-       * <tt>reinit()</tt> will have to
-       * give the vector the correct
-       * size and distribution among
-       * processes in case of an MPI
-       * run.
+       * Default constructor that generates an empty (zero size) vector. The
+       * function <tt>reinit()</tt> will have to give the vector the correct
+       * size and distribution among processes in case of an MPI run.
        */
       Vector ();
 
       /**
-       * Copy constructor using the
-       * given vector.
+       * Copy constructor using the given vector.
        */
       Vector (const Vector &V);
 
@@ -239,156 +244,88 @@ namespace TrilinosWrappers
       ~Vector ();
 
       /**
-       * Reinit functionality. This
-       * function sets the calling vector
-       * to the dimension and the parallel
-       * distribution of the input vector,
-       * but does not copy the elements in
-       * <tt>v</tt>. If <tt>fast</tt> is
-       * not <tt>true</tt>, the elements in
-       * the vector are initialized with
-       * zero, otherwise the content will
-       * be left unchanged and the user has
-       * to set all elements.
+       * Reinit functionality. This function sets the calling vector to the
+       * dimension and the parallel distribution of the input vector, but does
+       * not copy the elements in <tt>v</tt>. If <tt>fast</tt> is not
+       * <tt>true</tt>, the elements in the vector are initialized with zero,
+       * otherwise the content will be left unchanged and the user has to set
+       * all elements.
        *
-       * This function has a third  argument,
-       * <tt>allow_different_maps</tt>,
-       * that allows for an exchange of
-       * data between two equal-sized
-       * vectors (but being distributed
-       * differently among the
-       * processors). A trivial application
-       * of this function is to generate a
-       * replication of a whole vector on
-       * each machine, when the calling
-       * vector is built according to the
-       * localized vector class
-       * TrilinosWrappers::Vector, and
-       * <tt>v</tt> is a distributed
-       * vector. In this case, the variable
-       * <tt>fast</tt> needs to be set to
-       * <tt>false</tt>, since it does not
-       * make sense to exchange data
-       * between differently parallelized
-       * vectors without touching the
-       * elements.
+       * This function has a third argument, <tt>allow_different_maps</tt>,
+       * that allows for an exchange of data between two equal-sized vectors
+       * (but being distributed differently among the processors). A trivial
+       * application of this function is to generate a replication of a whole
+       * vector on each machine, when the calling vector is built according to
+       * the localized vector class TrilinosWrappers::Vector, and <tt>v</tt>
+       * is a distributed vector. In this case, the variable <tt>fast</tt>
+       * needs to be set to <tt>false</tt>, since it does not make sense to
+       * exchange data between differently parallelized vectors without
+       * touching the elements.
        */
       void reinit (const VectorBase &v,
                    const bool        fast = false,
                    const bool        allow_different_maps = false);
 
-      void reinit (const BlockVector &v,
-                   const bool         import_data = false);
-
       /**
-       * Reinit function. Creates a vector
-       * using the IndexSet local as our
-       * own unknowns, add optional ghost
-       * values ghost.
+       * Create vector by merging components from a block vector.
        */
-      void reinit (const IndexSet &local,
-                   const IndexSet &ghost,
-                   const MPI_Comm &communicator = MPI_COMM_WORLD);
+      void reinit (const BlockVector &v,
+                   const bool         import_data = false);
 
       /**
-       * Set all components of the
-       * vector to the given number @p
-       * s. Simply pass this down to
-       * the base class, but we still
-       * need to declare this function
-       * to make the example given in
-       * the discussion about making
-       * the constructor explicit work.
+       * Set all components of the vector to the given number @p s. Simply
+       * pass this down to the base class, but we still need to declare this
+       * function to make the example given in the discussion about making the
+       * constructor explicit work.
        */
       Vector &operator = (const TrilinosScalar s);
 
       /**
-       * Copy the given vector. Resize
-       * the present vector if
-       * necessary. In this case, also
-       * the Epetra_Map that designs
-       * the parallel partitioning is
-       * taken from the input vector.
+       * Copy the given vector. Resize the present vector if necessary. In
+       * this case, also the Epetra_Map that designs the parallel partitioning
+       * is taken from the input vector.
        */
       Vector &
       operator = (const Vector &V);
 
       /**
-       * Copy operator from a given
-       * localized vector (present on
-       * all processes) in
-       * TrilinosWrappers format to the
-       * current distributed
-       * vector. This function assumes
-       * that the calling vector (left
-       * hand object) already is of the
-       * same size as the right hand
-       * side vector. Otherwise, an
-       * exception will be thrown.
+       * Copy operator from a given localized vector (present on all
+       * processes) in TrilinosWrappers format to the current distributed
+       * vector. This function assumes that the calling vector (left hand
+       * object) already is of the same size as the right hand side
+       * vector. Otherwise, an exception will be thrown.
        */
       Vector &
       operator = (const ::dealii::TrilinosWrappers::Vector &V);
 
       /**
-       * Another copy function. This
-       * one takes a deal.II vector and
-       * copies it into a
-       * TrilinosWrapper vector. Note
-       * that since we do not provide
-       * any Epetra_map that tells
-       * about the partitioning of the
-       * vector among the MPI
-       * processes, the size of the
-       * TrilinosWrapper vector has to
-       * be the same as the size of the
-       * input vector. In order to
-       * change the map, use the
-       * reinit(const Epetra_Map
-       * &input_map) function.
+       * Another copy function. This one takes a deal.II vector and copies it
+       * into a TrilinosWrapper vector. Note that since we do not provide any
+       * Epetra_map that tells about the partitioning of the vector among the
+       * MPI processes, the size of the TrilinosWrapper vector has to be the
+       * same as the size of the input vector. In order to change the map, use
+       * the reinit(const Epetra_Map &input_map) function.
        */
       template <typename Number>
       Vector &
       operator = (const ::dealii::Vector<Number> &v);
 
       /**
-       * This reinit function is
-       * meant to be used for
-       * parallel calculations where
-       * some non-local data has to
-       * be used. The typical
-       * situation where one needs
-       * this function is the call of
-       * the
-       * FEValues<dim>::get_function_values
-       * function (or of some
-       * derivatives) in
-       * parallel. Since it is
-       * usually faster to retrieve
-       * the data in advance, this
-       * function can be called
-       * before the assembly forks
-       * out to the different
-       * processors. What this
-       * function does is the
-       * following: It takes the
-       * information in the columns
-       * of the given matrix and
-       * looks which data couples
-       * between the different
-       * processors. That data is
-       * then queried from the input
-       * vector. Note that you should
-       * not write to the resulting
-       * vector any more, since the
-       * some data can be stored
-       * several times on different
-       * processors, leading to
-       * unpredictable results. In
-       * particular, such a vector
-       * cannot be used for
-       * matrix-vector products as
-       * for example done during the
-       * solution of linear systems.
+       * This reinit function is meant to be used for parallel calculations
+       * where some non-local data has to be used. The typical situation where
+       * one needs this function is the call of the
+       * FEValues<dim>::get_function_values function (or of some derivatives)
+       * in parallel. Since it is usually faster to retrieve the data in
+       * advance, this function can be called before the assembly forks out to
+       * the different processors. What this function does is the following:
+       * It takes the information in the columns of the given matrix and looks
+       * which data couples between the different processors. That data is
+       * then queried from the input vector. Note that you should not write to
+       * the resulting vector any more, since the some data can be stored
+       * several times on different processors, leading to unpredictable
+       * results. In particular, such a vector cannot be used for
+       * matrix-vector products as for example done during the solution of
+       * linear systems.
        */
       void import_nonlocal_data_for_fe
       (const dealii::TrilinosWrappers::SparseMatrix &matrix,
@@ -399,35 +336,24 @@ namespace TrilinosWrappers
        */
 //@{
       /**
-       * This constructor takes an
-       * Epetra_Map that already knows
-       * how to distribute the
-       * individual components among
-       * the MPI processors. Since it
-       * also includes information
-       * about the size of the vector,
-       * this is all we need to
-       * generate a parallel vector.
+       * This constructor takes an Epetra_Map that already knows how to
+       * distribute the individual components among the MPI processors. Since
+       * it also includes information about the size of the vector, this is
+       * all we need to generate a parallel vector.
        */
       explicit Vector (const Epetra_Map &parallel_partitioning);
 
       /**
-       * Copy constructor from the
-       * TrilinosWrappers vector
-       * class. Since a vector of this
-       * class does not necessarily
-       * need to be distributed among
-       * processes, the user needs to
-       * supply us with an Epetra_Map
-       * that sets the partitioning
-       * details.
+       * Copy constructor from the TrilinosWrappers vector class. Since a
+       * vector of this class does not necessarily need to be distributed
+       * among processes, the user needs to supply us with an Epetra_Map that
+       * sets the partitioning details.
        */
       Vector (const Epetra_Map &parallel_partitioning,
               const VectorBase &v);
 
       /**
-       * Reinitialize from a deal.II
-       * vector. The Epetra_Map specifies the
+       * Reinitialize from a deal.II vector. The Epetra_Map specifies the
        * %parallel partitioning.
        */
       template <typename number>
@@ -435,20 +361,15 @@ namespace TrilinosWrappers
                    const dealii::Vector<number> &v);
 
       /**
-       * Reinit functionality. This
-       * function destroys the old
-       * vector content and generates a
-       * new one based on the input
-       * map.
+       * Reinit functionality. This function destroys the old vector content
+       * and generates a new one based on the input map.
        */
       void reinit (const Epetra_Map &parallel_partitioning,
                    const bool        fast = false);
 
       /**
-       * Copy-constructor from deal.II
-       * vectors. Sets the dimension to that
-       * of the given vector, and copies all
-       * elements.
+       * Copy-constructor from deal.II vectors. Sets the dimension to that of
+       * the given vector, and copies all elements.
        */
       template <typename Number>
       Vector (const Epetra_Map             &parallel_partitioning,
@@ -459,14 +380,10 @@ namespace TrilinosWrappers
        */
 //@{
       /**
-       * This constructor takes an IndexSet
-       * that defines how to distribute the
-       * individual components among the
-       * MPI processors. Since it also
-       * includes information about the
-       * size of the vector, this is all we
-       * need to generate a %parallel
-       * vector.
+       * This constructor takes an IndexSet that defines how to distribute the
+       * individual components among the MPI processors. Since it also
+       * includes information about the size of the vector, this is all we
+       * need to generate a %parallel vector.
        */
       explicit Vector (const IndexSet &parallel_partitioning,
                        const MPI_Comm &communicator = MPI_COMM_WORLD);
@@ -479,25 +396,18 @@ namespace TrilinosWrappers
               const MPI_Comm &communicator = MPI_COMM_WORLD);
 
       /**
-       * Copy constructor from the
-       * TrilinosWrappers vector
-       * class. Since a vector of this
-       * class does not necessarily need to
-       * be distributed among processes,
-       * the user needs to supply us with
-       * an IndexSet and an MPI
-       * communicator that set the
-       * partitioning details.
+       * Copy constructor from the TrilinosWrappers vector class. Since a
+       * vector of this class does not necessarily need to be distributed
+       * among processes, the user needs to supply us with an IndexSet and an
+       * MPI communicator that set the partitioning details.
        */
       Vector (const IndexSet   &parallel_partitioning,
               const VectorBase &v,
               const MPI_Comm   &communicator = MPI_COMM_WORLD);
 
       /**
-       * Copy-constructor from deal.II
-       * vectors. Sets the dimension to
-       * that of the given vector, and
-       * copies all the elements.
+       * Copy-constructor from deal.II vectors. Sets the dimension to that of
+       * the given vector, and copies all the elements.
        */
       template <typename Number>
       Vector (const IndexSet               &parallel_partitioning,
@@ -505,17 +415,37 @@ namespace TrilinosWrappers
               const MPI_Comm               &communicator = MPI_COMM_WORLD);
 
       /**
-       * Reinit functionality. This function
-       * destroys the old vector content and
-       * generates a new one based on the
-       * input partitioning.  The flag
-       * <tt>fast</tt> determines whether the
-       * vector should be filled with zero
-       * (false) or left untouched (true).
+       * Reinit functionality. This function destroys the old vector content
+       * and generates a new one based on the input partitioning.  The flag
+       * <tt>fast</tt> determines whether the vector should be filled with
+       * zero (false) or left untouched (true).
        */
       void reinit (const IndexSet &parallel_partitioning,
                    const MPI_Comm &communicator = MPI_COMM_WORLD,
                    const bool      fast = false);
+
+      /**
+       * Reinit functionality. This function destroys the old vector content
+       * and generates a new one based on the input partitioning. In addition
+       * to just specifying one index set as in all the other methods above,
+       * this method allows to supply an additional set of ghost
+       * entries. There are two different versions of a vector that can be
+       * created. If the flag @p vector_writable is set to @p false, the
+       * vector only allows read access to the joint set of @p
+       * parallel_partitioning and @p ghost_entries. The effect of the reinit
+       * method is then equivalent to calling the other reinit method with an
+       * index set containing both the locally owned entries and the ghost
+       * entries.
+       *
+       * If the flag @p vector_writable is set to true, this creates an
+       * alternative storage scheme for ghost elements that allows multiple
+       * threads to write into the vector (for the other reinit methods, only
+       * one thread is allowed to write into the ghost entries at a time).
+       */
+      void reinit (const IndexSet &locally_owned_entries,
+                   const IndexSet &ghost_entries,
+                   const MPI_Comm &communicator = MPI_COMM_WORLD,
+                   const bool      vector_writable = false);
 //@}
     };
 
@@ -639,124 +569,86 @@ namespace TrilinosWrappers
     typedef dealii::types::global_dof_index size_type;
 
     /**
-     * A variable that indicates whether this vector
-     * supports distributed data storage. If true, then
-     * this vector also needs an appropriate compress()
-     * function that allows communicating recent set or
-     * add operations to individual elements to be communicated
-     * to other processors.
+     * A variable that indicates whether this vector supports distributed data
+     * storage. If true, then this vector also needs an appropriate compress()
+     * function that allows communicating recent set or add operations to
+     * individual elements to be communicated to other processors.
      *
-     * For the current class, the variable equals
-     * false, since it does not support parallel data storage.
-     * If you do need parallel data storage, use
-     * TrilinosWrappers::MPI::Vector.
+     * For the current class, the variable equals false, since it does not
+     * support parallel data storage.  If you do need parallel data storage,
+     * use TrilinosWrappers::MPI::Vector.
      */
     static const bool supports_distributed_data = false;
 
     /**
-     * Default constructor that
-     * generates an empty (zero size)
-     * vector. The function
-     * <tt>reinit()</tt> will have to
-     * give the vector the correct
+     * Default constructor that generates an empty (zero size) vector. The
+     * function <tt>reinit()</tt> will have to give the vector the correct
      * size.
      */
     Vector ();
 
     /**
-     * This constructor takes as
-     * input the number of elements
-     * in the vector.
+     * This constructor takes as input the number of elements in the vector.
      */
     explicit Vector (const size_type n);
 
     /**
-     * This constructor takes as
-     * input the number of elements
-     * in the vector. If the map is
-     * not localized, i.e., if there
-     * are some elements that are not
-     * present on all processes, only
-     * the global size of the map
-     * will be taken and a localized
-     * map will be generated
-     * internally.
+     * This constructor takes as input the number of elements in the
+     * vector. If the map is not localized, i.e., if there are some elements
+     * that are not present on all processes, only the global size of the map
+     * will be taken and a localized map will be generated internally.
      */
     explicit Vector (const Epetra_Map &partitioning);
 
     /**
-     * This constructor takes as input
-     * the number of elements in the
-     * vector. If the index set is not
-     * localized, i.e., if there are some
-     * elements that are not present on
-     * all processes, only the global
-     * size of the index set will be
-     * taken and a localized version will
-     * be generated internally.
+     * This constructor takes as input the number of elements in the
+     * vector. If the index set is not localized, i.e., if there are some
+     * elements that are not present on all processes, only the global size of
+     * the index set will be taken and a localized version will be generated
+     * internally.
      */
     explicit Vector (const IndexSet &partitioning,
                      const MPI_Comm &communicator = MPI_COMM_WORLD);
 
     /**
-     * This constructor takes a
-     * (possibly parallel) Trilinos
-     * Vector and generates a
-     * localized version of the whole
-     * content on each processor.
+     * This constructor takes a (possibly parallel) Trilinos Vector and
+     * generates a localized version of the whole content on each processor.
      */
     explicit Vector (const VectorBase &V);
 
     /**
-     * Copy-constructor from deal.II
-     * vectors. Sets the dimension to that
-     * of the given vector, and copies all
-     * elements.
+     * Copy-constructor from deal.II vectors. Sets the dimension to that of
+     * the given vector, and copies all elements.
      */
     template <typename Number>
     explicit Vector (const dealii::Vector<Number> &v);
 
     /**
-     * Reinit function that resizes
-     * the vector to the size
-     * specified by <tt>n</tt>.
+     * Reinit function that resizes the vector to the size specified by
+     * <tt>n</tt>.
      */
     void reinit (const size_type n,
                  const bool      fast = false);
 
     /**
-     * Initialization with an
-     * Epetra_Map. Similar to the call in
-     * the other class MPI::Vector, with
-     * the difference that now a copy on
-     * all processes is generated. This
-     * initialization function is
-     * appropriate when the data in the
-     * localized vector should be
-     * imported from a distributed vector
-     * that has been initialized with the
-     * same communicator. The variable
-     * <tt>fast</tt> determines whether
-     * the vector should be filled with
-     * zero or left untouched.
+     * Initialization with an Epetra_Map. Similar to the call in the other
+     * class MPI::Vector, with the difference that now a copy on all processes
+     * is generated. This initialization function is appropriate when the data
+     * in the localized vector should be imported from a distributed vector
+     * that has been initialized with the same communicator. The variable
+     * <tt>fast</tt> determines whether the vector should be filled with zero
+     * or left untouched.
      */
     void reinit (const Epetra_Map &input_map,
                  const bool        fast = false);
 
     /**
-     * Initialization with an
-     * IndexSet. Similar to the call in the
-     * other class MPI::Vector, with the
-     * difference that now a copy on all
-     * processes is generated. This
-     * initialization function is
-     * appropriate in case the data in the
-     * localized vector should be imported
-     * from a distributed vector that has
-     * been initialized with the same
-     * communicator. The variable
-     * <tt>fast</tt> determines whether the
-     * vector should be filled with zero
+     * Initialization with an IndexSet. Similar to the call in the other class
+     * MPI::Vector, with the difference that now a copy on all processes is
+     * generated. This initialization function is appropriate in case the data
+     * in the localized vector should be imported from a distributed vector
+     * that has been initialized with the same communicator. The variable
+     * <tt>fast</tt> determines whether the vector should be filled with zero
      * (false) or left untouched (true).
      */
     void reinit (const IndexSet   &input_map,
@@ -764,66 +656,52 @@ namespace TrilinosWrappers
                  const bool        fast = false);
 
     /**
-     * Reinit function. Takes the
-     * information of a Vector and copies
-     * everything to the calling vector,
-     * now also allowing different maps.
+     * Reinit function. Takes the information of a Vector and copies
+     * everything to the calling vector, now also allowing different maps.
      */
     void reinit (const VectorBase &V,
                  const bool        fast = false,
                  const bool        allow_different_maps = false);
 
     /**
-     * Set all components of the
-     * vector to the given number @p
-     * s. Simply pass this down to
-     * the base class, but we still
-     * need to declare this function
-     * to make the example given in
-     * the discussion about making
-     * the constructor explicit work.
+     * Set all components of the vector to the given number @p s. Simply pass
+     * this down to the base class, but we still need to declare this function
+     * to make the example given in the discussion about making the
+     * constructor explicit work.
      */
     Vector &operator = (const TrilinosScalar s);
 
     /**
-     * Sets the left hand argument to
-     * the (parallel) Trilinos
-     * Vector. Equivalent to the @p
-     * reinit function.
+     * Sets the left hand argument to the (parallel) Trilinos
+     * Vector. Equivalent to the @p reinit function.
      */
     Vector &
     operator = (const MPI::Vector &V);
 
     /**
-     * Sets the left hand argument to
-     * the deal.II vector.
+     * Sets the left hand argument to the deal.II vector.
      */
     template <typename Number>
     Vector &
     operator = (const ::dealii::Vector<Number> &V);
 
     /**
-     * Copy operator. Copies both the
-     * dimension and the content in
-     * the right hand argument.
+     * Copy operator. Copies both the dimension and the content in the right
+     * hand argument.
      */
     Vector &
     operator = (const Vector &V);
 
     /**
-     * This function does nothing but is
-     * there for compatibility with the
-     * @p PETScWrappers::Vector class.
+     * This function does nothing but is there for compatibility with the @p
+     * PETScWrappers::Vector class.
      *
-     * For the PETSc vector wrapper class,
-     * this function updates the ghost
-     * values of the PETSc vector. This
-     * is necessary after any modification
+     * For the PETSc vector wrapper class, this function updates the ghost
+     * values of the PETSc vector. This is necessary after any modification
      * before reading ghost values.
      *
-     * However, for the implementation of
-     * this class, it is immaterial and thus
-     * an empty function.
+     * However, for the implementation of this class, it is immaterial and
+     * thus an empty function.
      */
     void update_ghost_values () const;
   };
index 66995f19cdc04203431d936c21129d4340f30e0e..39b9b79991326845a0c8fe756c3cc79273eb07b1 100644 (file)
@@ -75,29 +75,20 @@ namespace TrilinosWrappers
     typedef dealii::types::global_dof_index size_type;
 
     /**
-     * This class implements a
-     * wrapper for accessing the
-     * Trilinos vector in the same
-     * way as we access deal.II
-     * objects: it is initialized
-     * with a vector and an element
-     * within it, and has a
-     * conversion operator to
-     * extract the scalar value of
-     * this element. It also has a
-     * variety of assignment
-     * operator for writing to this
-     * one element.  @ingroup
-     * TrilinosWrappers
+     * This class implements a wrapper for accessing the Trilinos vector in
+     * the same way as we access deal.II objects: it is initialized with a
+     * vector and an element within it, and has a conversion operator to
+     * extract the scalar value of this element. It also has a variety of
+     * assignment operator for writing to this one element.
+     *
+     * @ingroup TrilinosWrappers
      */
     class VectorReference
     {
     private:
       /**
-       * Constructor. It is made
-       * private so as to only allow
-       * the actual vector class to
-       * create it.
+       * Constructor. It is made private so as to only allow the actual vector
+       * class to create it.
        */
       VectorReference (VectorBase     &vector,
                        const size_type  index);
@@ -105,83 +96,58 @@ namespace TrilinosWrappers
     public:
 
       /**
-       * This looks like a copy
-       * operator, but does something
-       * different than usual. In
-       * particular, it does not copy
-       * the member variables of this
-       * reference. Rather, it
-       * handles the situation where
-       * we have two vectors @p v and
-       * @p w, and assign elements
-       * like in
-       * <tt>v(i)=w(i)</tt>. Here,
-       * both left and right hand
-       * side of the assignment have
-       * data type VectorReference,
-       * but what we really mean is
-       * to assign the vector
-       * elements represented by the
-       * two references. This
-       * operator implements this
-       * operation. Note also that
-       * this allows us to make the
-       * assignment operator const.
+       * This looks like a copy operator, but does something different than
+       * usual. In particular, it does not copy the member variables of this
+       * reference. Rather, it handles the situation where we have two vectors
+       * @p v and @p w, and assign elements like in <tt>v(i)=w(i)</tt>. Here,
+       * both left and right hand side of the assignment have data type
+       * VectorReference, but what we really mean is to assign the vector
+       * elements represented by the two references. This operator implements
+       * this operation. Note also that this allows us to make the assignment
+       * operator const.
        */
       const VectorReference &
       operator = (const VectorReference &r) const;
 
       /**
-       * Same as above but for non-const
-       * reference objects.
+       * Same as above but for non-const reference objects.
        */
       const VectorReference &
       operator = (const VectorReference &r);
 
       /**
-       * Set the referenced element of the
-       * vector to <tt>s</tt>.
+       * Set the referenced element of the vector to <tt>s</tt>.
        */
       const VectorReference &
       operator = (const TrilinosScalar &s) const;
 
       /**
-       * Add <tt>s</tt> to the
-       * referenced element of the
-       * vector->
+       * Add <tt>s</tt> to the referenced element of the vector->
        */
       const VectorReference &
       operator += (const TrilinosScalar &s) const;
 
       /**
-       * Subtract <tt>s</tt> from the
-       * referenced element of the
-       * vector->
+       * Subtract <tt>s</tt> from the referenced element of the vector->
        */
       const VectorReference &
       operator -= (const TrilinosScalar &s) const;
 
       /**
-       * Multiply the referenced
-       * element of the vector by
-       * <tt>s</tt>.
+       * Multiply the referenced element of the vector by <tt>s</tt>.
        */
       const VectorReference &
       operator *= (const TrilinosScalar &s) const;
 
       /**
-       * Divide the referenced
-       * element of the vector by
-       * <tt>s</tt>.
+       * Divide the referenced element of the vector by <tt>s</tt>.
        */
       const VectorReference &
       operator /= (const TrilinosScalar &s) const;
 
       /**
-       * Convert the reference to an
-       * actual value, i.e. return
-       * the value of the referenced
-       * element of the vector.
+       * Convert the reference to an actual value, i.e. return the value of
+       * the referenced element of the vector.
        */
       operator TrilinosScalar () const;
 
@@ -209,21 +175,17 @@ namespace TrilinosWrappers
 
     private:
       /**
-       * Point to the vector we are
-       * referencing.
+       * Point to the vector we are referencing.
        */
       VectorBase   &vector;
 
       /**
-       * Index of the referenced element
-       * of the vector.
+       * Index of the referenced element of the vector.
        */
       const size_type  index;
 
       /**
-       * Make the vector class a
-       * friend, so that it can
-       * create objects of the
+       * Make the vector class a friend, so that it can create objects of the
        * present type.
        */
       friend class ::dealii::TrilinosWrappers::VectorBase;
@@ -270,11 +232,8 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Declare some of the standard
-     * types used in all
-     * containers. These types
-     * parallel those in the
-     * <tt>C</tt> standard libraries
+     * Declare some of the standard types used in all containers. These types
+     * parallel those in the <tt>C</tt> standard libraries
      * <tt>vector<...></tt> class.
      */
     typedef TrilinosScalar                  value_type;
@@ -291,22 +250,15 @@ namespace TrilinosWrappers
     //@{
 
     /**
-     * Default constructor that
-     * generates an empty (zero size)
-     * vector. The function
-     * <tt>reinit()</tt> will have to
-     * give the vector the correct
-     * size and distribution among
-     * processes in case of an MPI
-     * run.
+     * Default constructor that generates an empty (zero size) vector. The
+     * function <tt>reinit()</tt> will have to give the vector the correct
+     * size and distribution among processes in case of an MPI run.
      */
     VectorBase ();
 
     /**
-     * Copy constructor. Sets the
-     * dimension to that of the given
-     * vector, and copies all the
-     * elements.
+     * Copy constructor. Sets the dimension to that of the given vector, and
+     * copies all the elements.
      */
     VectorBase (const VectorBase &v);
 
@@ -316,47 +268,30 @@ namespace TrilinosWrappers
     virtual ~VectorBase ();
 
     /**
-     * Release all memory and return
-     * to a state just like after
-     * having called the default
-     * constructor.
+     * Release all memory and return to a state just like after having called
+     * the default constructor.
      */
     void clear ();
 
     /**
-     * Reinit functionality, sets the
-     * dimension and possibly the
-     * parallel partitioning (Epetra_Map)
-     * of the calling vector to the
-     * settings of the input vector.
+     * Reinit functionality, sets the dimension and possibly the parallel
+     * partitioning (Epetra_Map) of the calling vector to the settings of the
+     * input vector.
      */
     void reinit (const VectorBase &v,
                  const bool        fast = false);
 
     /**
-     * Compress the underlying
-     * representation of the Trilinos
-     * object, i.e. flush the buffers
-     * of the vector object if it has
-     * any. This function is
-     * necessary after writing into a
-     * vector element-by-element and
-     * before anything else can be
-     * done on it.
+     * Compress the underlying representation of the Trilinos object,
+     * i.e. flush the buffers of the vector object if it has any. This
+     * function is necessary after writing into a vector element-by-element
+     * and before anything else can be done on it.
      *
-     * The (defaulted) argument can
-     * be used to specify the
-     * compress mode
-     * (<code>Add</code> or
-     * <code>Insert</code>) in case
-     * the vector has not been
-     * written to since the last
-     * time this function was
-     * called. The argument is
-     * ignored if the vector has
-     * been added or written to
-     * since the last time
-     * compress() was called.
+     * The (defaulted) argument can be used to specify the compress mode
+     * (<code>Add</code> or <code>Insert</code>) in case the vector has not
+     * been written to since the last time this function was called. The
+     * argument is ignored if the vector has been added or written to since
+     * the last time compress() was called.
      *
      * See @ref GlossCompress "Compressing distributed objects"
      * for more information.
@@ -375,136 +310,87 @@ namespace TrilinosWrappers
     void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED;
 
     /**
-     * Returns the state of the
-     * vector, i.e., whether
-     * compress() has already been
-     * called after an operation
-     * requiring data exchange.
+     * Returns the state of the vector, i.e., whether compress() has already
+     * been called after an operation requiring data exchange.
      */
     bool is_compressed () const;
 
     /**
-     * Set all components of the
-     * vector to the given number @p
-     * s. Simply pass this down to
-     * the Trilinos Epetra object,
-     * but we still need to declare
-     * this function to make the
-     * example given in the
-     * discussion about making the
-     * constructor explicit work.
+     * Set all components of the vector to the given number @p s. Simply pass
+     * this down to the Trilinos Epetra object, but we still need to declare
+     * this function to make the example given in the discussion about making
+     * the constructor explicit work.
      *
-     * Since the semantics of
-     * assigning a scalar to a vector
-     * are not immediately clear,
-     * this operator should really
-     * only be used if you want to
-     * set the entire vector to
-     * zero. This allows the
-     * intuitive notation
-     * <tt>v=0</tt>. Assigning other
-     * values is deprecated and may
-     * be disallowed in the future.
+     * Since the semantics of assigning a scalar to a vector are not
+     * immediately clear, this operator should really only be used if you want
+     * to set the entire vector to zero. This allows the intuitive notation
+     * <tt>v=0</tt>. Assigning other values is deprecated and may be
+     * disallowed in the future.
      */
     VectorBase &
     operator = (const TrilinosScalar s);
 
     /**
-     * Copy function. This function takes
-     * a VectorBase vector and copies all
-     * the elements. The target vector
-     * will have the same parallel
-     * distribution as the calling
-     * vector.
+     * Copy function. This function takes a VectorBase vector and copies all
+     * the elements. The target vector will have the same parallel
+     * distribution as the calling vector.
      */
     VectorBase &
     operator = (const VectorBase &v);
 
     /**
-     * Another copy function. This
-     * one takes a deal.II vector and
-     * copies it into a
-     * TrilinosWrapper vector. Note
-     * that since we do not provide
-     * any Epetra_map that tells
-     * about the partitioning of the
-     * vector among the MPI
-     * processes, the size of the
-     * TrilinosWrapper vector has to
-     * be the same as the size of the
-     * input vector. In order to
-     * change the map, use the
-     * reinit(const Epetra_Map
-     * &input_map) function.
+     * Another copy function. This one takes a deal.II vector and copies it
+     * into a TrilinosWrapper vector. Note that since we do not provide any
+     * Epetra_map that tells about the partitioning of the vector among the
+     * MPI processes, the size of the TrilinosWrapper vector has to be the
+     * same as the size of the input vector. In order to change the map, use
+     * the reinit(const Epetra_Map &input_map) function.
      */
     template <typename Number>
     VectorBase &
     operator = (const ::dealii::Vector<Number> &v);
 
     /**
-     * Test for equality. This
-     * function assumes that the
-     * present vector and the one to
-     * compare with have the same
-     * size already, since comparing
-     * vectors of different sizes
-     * makes not much sense anyway.
+     * Test for equality. This function assumes that the present vector and
+     * the one to compare with have the same size already, since comparing
+     * vectors of different sizes makes not much sense anyway.
      */
     bool operator == (const VectorBase &v) const;
 
     /**
-     * Test for inequality. This
-     * function assumes that the
-     * present vector and the one to
-     * compare with have the same
-     * size already, since comparing
-     * vectors of different sizes
-     * makes not much sense anyway.
+     * Test for inequality. This function assumes that the present vector and
+     * the one to compare with have the same size already, since comparing
+     * vectors of different sizes makes not much sense anyway.
      */
     bool operator != (const VectorBase &v) const;
 
     /**
-     * Return the global dimension of
-     * the vector.
+     * Return the global dimension of the vector.
      */
     size_type size () const;
 
     /**
-     * Return the local dimension of
-     * the vector, i.e. the number of
-     * elements stored on the present
-     * MPI process. For sequential
-     * vectors, this number is the
-     * same as size(), but for
-     * parallel vectors it may be
-     * smaller.
+     * Return the local dimension of the vector, i.e. the number of elements
+     * stored on the present MPI process. For sequential vectors, this number
+     * is the same as size(), but for parallel vectors it may be smaller.
      *
-     * To figure out which elements
-     * exactly are stored locally,
-     * use local_range().
+     * To figure out which elements exactly are stored locally, use
+     * local_range().
      *
-     * If the vector contains ghost
-     * elements, they are included in
-     * this number.
+     * If the vector contains ghost elements, they are included in this
+     * number.
      */
     size_type local_size () const;
 
     /**
-     * Return a pair of indices
-     * indicating which elements of
-     * this vector are stored
-     * locally. The first number is
-     * the index of the first element
-     * stored, the second the index
-     * of the one past the last one
-     * that is stored locally. If
-     * this is a sequential vector,
-     * then the result will be the
-     * pair <code>(0,N)</code>, otherwise it will
-     * be a pair <code>(i,i+n)</code>, where
-     * <code>n=local_size()</code> and <code>i</code> is the first
-     * element of the vector stored on this processor, corresponding
-     * to the half open interval $[i,i+n)$
+     * Return a pair of indices indicating which elements of this vector are
+     * stored locally. The first number is the index of the first element
+     * stored, the second the index of the one past the last one that is
+     * stored locally. If this is a sequential vector, then the result will be
+     * the pair <code>(0,N)</code>, otherwise it will be a pair
+     * <code>(i,i+n)</code>, where <code>n=local_size()</code> and
+     * <code>i</code> is the first element of the vector stored on this
+     * processor, corresponding to the half open interval $[i,i+n)$
      *
      * @note The description above is true most of the time, but
      * not always. In particular, Trilinos vectors need not store
@@ -522,9 +408,8 @@ namespace TrilinosWrappers
     std::pair<size_type, size_type> local_range () const;
 
     /**
-     * Return whether @p index is in
-     * the local range or not, see
-     * also local_range().
+     * Return whether @p index is in the local range or not, see also
+     * local_range().
      *
      * @note The same limitation for the applicability of this
      * function applies as listed in the documentation of local_range().
@@ -548,88 +433,65 @@ namespace TrilinosWrappers
     IndexSet locally_owned_elements () const;
 
     /**
-     * Return if the vector contains ghost
-     * elements. This answer is true if there
-     * are ghost elements on at least one
-     * process.
+     * Return if the vector contains ghost elements. This answer is true if
+     * there are ghost elements on at least one process.
      */
     bool has_ghost_elements() const;
 
     /**
-     * Return the scalar (inner)
-     * product of two vectors. The
-     * vectors must have the same
-     * size.
+     * Return the scalar (inner) product of two vectors. The vectors must have
+     * the same size.
      */
     TrilinosScalar operator * (const VectorBase &vec) const;
 
     /**
-     * Return square of the
-     * $l_2$-norm.
+     * Return square of the $l_2$-norm.
      */
     real_type norm_sqr () const;
 
     /**
-     * Mean value of the elements of
-     * this vector.
+     * Mean value of the elements of this vector.
      */
     TrilinosScalar mean_value () const;
 
     /**
-     * Compute the minimal value of
-     * the elements of this vector.
+     * Compute the minimal value of the elements of this vector.
      */
     TrilinosScalar minimal_value () const;
 
     /**
-     * $l_1$-norm of the vector.  The
-     * sum of the absolute values.
+     * $l_1$-norm of the vector.  The sum of the absolute values.
      */
     real_type l1_norm () const;
 
     /**
-     * $l_2$-norm of the vector.  The
-     * square root of the sum of the
-     * squares of the elements.
+     * $l_2$-norm of the vector.  The square root of the sum of the squares of
+     * the elements.
      */
     real_type l2_norm () const;
 
     /**
-     * $l_p$-norm of the vector. The
-     * <i>p</i>th root of the sum of
-     * the <i>p</i>th powers of the
-     * absolute values of the
-     * elements.
+     * $l_p$-norm of the vector. The <i>p</i>th root of the sum of the
+     * <i>p</i>th powers of the absolute values of the elements.
      */
     real_type lp_norm (const TrilinosScalar p) const;
 
     /**
-     * Maximum absolute value of the
-     * elements.
+     * Maximum absolute value of the elements.
      */
     real_type linfty_norm () const;
 
     /**
-     * Return whether the vector
-     * contains only elements with
-     * value zero. This function is
-     * mainly for internal
-     * consistency checks and should
-     * seldom be used when not in
-     * debug mode since it uses quite
-     * some time.
+     * Return whether the vector contains only elements with value zero. This
+     * function is mainly for internal consistency checks and should seldom be
+     * used when not in debug mode since it uses quite some time.
      */
     bool all_zero () const;
 
     /**
-     * Return @p true if the vector
-     * has no negative entries,
-     * i.e. all entries are zero or
-     * positive. This function is
-     * used, for example, to check
-     * whether refinement indicators
-     * are really all positive (or
-     * zero).
+     * Return @p true if the vector has no negative entries, i.e. all entries
+     * are zero or positive. This function is used, for example, to check
+     * whether refinement indicators are really all positive (or zero).
      */
     bool is_non_negative () const;
     //@}
@@ -641,23 +503,20 @@ namespace TrilinosWrappers
     //@{
 
     /**
-     * Provide access to a given
-     * element, both read and write.
+     * Provide access to a given element, both read and write.
      */
     reference
     operator () (const size_type index);
 
     /**
-     * Provide read-only access to an
-     * element. This is equivalent to
-     * the <code>el()</code> command.
+     * Provide read-only access to an element. This is equivalent to the
+     * <code>el()</code> command.
      */
     TrilinosScalar
     operator () (const size_type index) const;
 
     /**
-     * Provide access to a given
-     * element, both read and write.
+     * Provide access to a given element, both read and write.
      *
      * Exactly the same as operator().
      */
@@ -665,9 +524,8 @@ namespace TrilinosWrappers
     operator [] (const size_type index);
 
     /**
-     * Provide read-only access to an
-     * element. This is equivalent to
-     * the <code>el()</code> command.
+     * Provide read-only access to an element. This is equivalent to the
+     * <code>el()</code> command.
      *
      * Exactly the same as operator().
      */
@@ -675,21 +533,17 @@ namespace TrilinosWrappers
     operator [] (const size_type index) const;
 
     /**
-     * A collective get operation: instead
-     * of getting individual elements of a
-     * vector, this function allows to get
-     * a whole set of elements at once. The
-     * indices of the elements to be read
-     * are stated in the first argument,
-     * the corresponding values are returned in the
-     * second.
+     * A collective get operation: instead of getting individual elements of a
+     * vector, this function allows to get a whole set of elements at
+     * once. The indices of the elements to be read are stated in the first
+     * argument, the corresponding values are returned in the second.
      */
     void extract_subvector_to (const std::vector<size_type> &indices,
                                std::vector<TrilinosScalar> &values) const;
 
     /**
-     * Just as the above, but with pointers.
-     * Useful in minimizing copying of data around.
+     * Just as the above, but with pointers.  Useful in minimizing copying of
+     * data around.
      */
     template <typename ForwardIterator, typename OutputIterator>
     void extract_subvector_to (ForwardIterator          indices_begin,
@@ -697,22 +551,20 @@ namespace TrilinosWrappers
                                OutputIterator           values_begin) const;
 
     /**
-     * Return the value of the vector
-     * entry <i>i</i>. Note that this
-     * function does only work
-     * properly when we request a
-     * data stored on the local
-     * processor. The function will
-     * throw an exception in case the
-     * elements sits on another
-     * process.
+     * Return the value of the vector entry <i>i</i>. Note that this function
+     * does only work properly when we request a data stored on the local
+     * processor. The function will throw an exception in case the elements
+     * sits on another process.
      */
     TrilinosScalar el (const size_type index) const;
 
     /**
-     * Make the Vector class a bit like the <tt>vector<></tt> class of
-     * the C++ standard library by returning iterators to the start and end
-     * of the locally owned elements of this vector. The ordering of local elements corresponds to the one given
+     * Make the Vector class a bit like the <tt>vector<></tt> class of the C++
+     * standard library by returning iterators to the start and end of the
+     * locally owned elements of this vector. The ordering of local elements
+     * corresponds to the one given by the global indices in case the vector
+     * is constructed from an IndexSet or other methods in deal.II (note that
+     * an Epetra_Map can contain elements in arbitrary orders, though).
      *
      * It holds that end() - begin() == local_size().
      */
@@ -731,31 +583,23 @@ namespace TrilinosWrappers
     iterator end ();
 
     /**
-     * Return a constant iterator pointing to the element past the end of
-     * the array of the locally owned entries.
+     * Return a constant iterator pointing to the element past the end of the
+     * array of the locally owned entries.
      */
     const_iterator end () const;
 
     /**
-     * A collective set operation:
-     * instead of setting individual
-     * elements of a vector, this
-     * function allows to set a whole
-     * set of elements at once. The
-     * indices of the elements to be
-     * set are stated in the first
-     * argument, the corresponding
-     * values in the second.
+     * A collective set operation: instead of setting individual elements of a
+     * vector, this function allows to set a whole set of elements at
+     * once. The indices of the elements to be set are stated in the first
+     * argument, the corresponding values in the second.
      */
     void set (const std::vector<size_type>    &indices,
               const std::vector<TrilinosScalar>  &values);
 
     /**
-     * This is a second collective
-     * set operation. As a
-     * difference, this function
-     * takes a deal.II vector of
-     * values.
+     * This is a second collective set operation. As a difference, this
+     * function takes a deal.II vector of values.
      */
     void set (const std::vector<size_type>        &indices,
               const ::dealii::Vector<TrilinosScalar> &values);
@@ -768,109 +612,80 @@ namespace TrilinosWrappers
     //@{
 
     /**
-     * This collective set operation
-     * is of lower level and can
-     * handle anything else &mdash;
-     * the only thing you have to
-     * provide is an address where
-     * all the indices are stored and
-     * the number of elements to be
-     * set.
+     * This collective set operation is of lower level and can handle anything
+     * else &mdash; the only thing you have to provide is an address where all
+     * the indices are stored and the number of elements to be set.
      */
     void set (const size_type       n_elements,
               const size_type      *indices,
               const TrilinosScalar *values);
 
     /**
-     * A collective add operation:
-     * This funnction adds a whole
-     * set of values stored in @p
-     * values to the vector
-     * components specified by @p
-     * indices.
+     * A collective add operation: This funnction adds a whole set of values
+     * stored in @p values to the vector components specified by @p indices.
      */
     void add (const std::vector<size_type>      &indices,
               const std::vector<TrilinosScalar> &values);
 
     /**
-     * This is a second collective
-     * add operation. As a
-     * difference, this function
-     * takes a deal.II vector of
-     * values.
+     * This is a second collective add operation. As a difference, this
+     * function takes a deal.II vector of values.
      */
     void add (const std::vector<size_type>           &indices,
               const ::dealii::Vector<TrilinosScalar> &values);
 
     /**
-     * Take an address where
-     * <tt>n_elements</tt> are stored
-     * contiguously and add them into
-     * the vector. Handles all cases
-     * which are not covered by the
-     * other two <tt>add()</tt>
-     * functions above.
+     * Take an address where <tt>n_elements</tt> are stored contiguously and
+     * add them into the vector. Handles all cases which are not covered by
+     * the other two <tt>add()</tt> functions above.
      */
     void add (const size_type       n_elements,
               const size_type      *indices,
               const TrilinosScalar *values);
 
     /**
-     * Multiply the entire vector by
-     * a fixed factor.
+     * Multiply the entire vector by a fixed factor.
      */
     VectorBase &operator *= (const TrilinosScalar factor);
 
     /**
-     * Divide the entire vector by a
-     * fixed factor.
+     * Divide the entire vector by a fixed factor.
      */
     VectorBase &operator /= (const TrilinosScalar factor);
 
     /**
-     * Add the given vector to the
-     * present one.
+     * Add the given vector to the present one.
      */
     VectorBase &operator += (const VectorBase &V);
 
     /**
-     * Subtract the given vector from
-     * the present one.
+     * Subtract the given vector from the present one.
      */
     VectorBase &operator -= (const VectorBase &V);
 
     /**
-     * Addition of @p s to all
-     * components. Note that @p s is
-     * a scalar and not a vector.
+     * Addition of @p s to all components. Note that @p s is a scalar and not
+     * a vector.
      */
     void add (const TrilinosScalar s);
 
     /**
-     * Simple vector addition, equal
-     * to the <tt>operator
-     * +=</tt>.
+     * Simple vector addition, equal to the <tt>operator +=</tt>.
      *
-     * Though, if the second argument
-     * <tt>allow_different_maps</tt>
-     * is set, then it is possible to
-     * add data from a different map.
+     * Though, if the second argument <tt>allow_different_maps</tt> is set,
+     * then it is possible to add data from a different map.
      */
     void add (const VectorBase &V,
               const bool        allow_different_maps = false);
 
     /**
-     * Simple addition of a multiple
-     * of a vector, i.e. <tt>*this =
-     * a*V</tt>.
+     * Simple addition of a multiple of a vector, i.e. <tt>*this = a*V</tt>.
      */
     void add (const TrilinosScalar  a,
               const VectorBase     &V);
 
     /**
-     * Multiple addition of scaled
-     * vectors, i.e. <tt>*this = a*V +
-     * b*W</tt>.
+     * Multiple addition of scaled vectors, i.e. <tt>*this = a*V + b*W</tt>.
      */
     void add (const TrilinosScalar  a,
               const VectorBase     &V,
@@ -878,17 +693,14 @@ namespace TrilinosWrappers
               const VectorBase     &W);
 
     /**
-     * Scaling and simple vector
-     * addition, i.e.  <tt>*this =
-     * s*(*this) + V</tt>.
+     * Scaling and simple vector addition, i.e.  <tt>*this = s*(*this) +
+     * V</tt>.
      */
     void sadd (const TrilinosScalar  s,
                const VectorBase     &V);
 
     /**
-     * Scaling and simple addition,
-     * i.e.  <tt>*this = s*(*this) +
-     * a*V</tt>.
+     * Scaling and simple addition, i.e.  <tt>*this = s*(*this) + a*V</tt>.
      */
     void sadd (const TrilinosScalar  s,
                const TrilinosScalar  a,
@@ -904,9 +716,8 @@ namespace TrilinosWrappers
                const VectorBase     &W);
 
     /**
-     * Scaling and multiple addition.
-     * <tt>*this = s*(*this) + a*V +
-     * b*W + c*X</tt>.
+     * Scaling and multiple addition.  <tt>*this = s*(*this) + a*V + b*W +
+     * c*X</tt>.
      */
     void sadd (const TrilinosScalar  s,
                const TrilinosScalar  a,
@@ -917,26 +728,20 @@ namespace TrilinosWrappers
                const VectorBase     &X);
 
     /**
-     * Scale each element of this
-     * vector by the corresponding
-     * element in the argument. This
-     * function is mostly meant to
-     * simulate multiplication (and
-     * immediate re-assignment) by a
-     * diagonal scaling matrix.
+     * Scale each element of this vector by the corresponding element in the
+     * argument. This function is mostly meant to simulate multiplication (and
+     * immediate re-assignment) by a diagonal scaling matrix.
      */
     void scale (const VectorBase &scaling_factors);
 
     /**
-     * Assignment <tt>*this =
-     * a*V</tt>.
+     * Assignment <tt>*this = a*V</tt>.
      */
     void equ (const TrilinosScalar  a,
               const VectorBase     &V);
 
     /**
-     * Assignment <tt>*this = a*V +
-     * b*W</tt>.
+     * Assignment <tt>*this = a*V + b*W</tt>.
      */
     void equ (const TrilinosScalar  a,
               const VectorBase     &V,
@@ -944,21 +749,14 @@ namespace TrilinosWrappers
               const VectorBase     &W);
 
     /**
-     * Compute the elementwise ratio
-     * of the two given vectors, that
-     * is let <tt>this[i] =
-     * a[i]/b[i]</tt>. This is useful
-     * for example if you want to
-     * compute the cellwise ratio of
-     * true to estimated error.
+     * Compute the elementwise ratio of the two given vectors, that is let
+     * <tt>this[i] = a[i]/b[i]</tt>. This is useful for example if you want to
+     * compute the cellwise ratio of true to estimated error.
      *
-     * This vector is appropriately
-     * scaled to hold the result.
+     * This vector is appropriately scaled to hold the result.
      *
-     * If any of the <tt>b[i]</tt> is
-     * zero, the result is
-     * undefined. No attempt is made
-     * to catch such situations.
+     * If any of the <tt>b[i]</tt> is zero, the result is undefined. No
+     * attempt is made to catch such situations.
      */
     void ratio (const VectorBase &a,
                 const VectorBase &b);
@@ -971,46 +769,34 @@ namespace TrilinosWrappers
     //@{
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos
-     * Epetra_MultiVector class.
+     * Return a const reference to the underlying Trilinos Epetra_MultiVector
+     * class.
      */
     const Epetra_MultiVector &trilinos_vector () const;
 
     /**
-     * Return a (modifyable) reference to
-     * the underlying Trilinos
+     * Return a (modifyable) reference to the underlying Trilinos
      * Epetra_FEVector class.
      */
     Epetra_FEVector &trilinos_vector ();
 
     /**
-     * Return a const reference to the
-     * underlying Trilinos Epetra_Map
-     * that sets the parallel
-     * partitioning of the vector.
+     * Return a const reference to the underlying Trilinos Epetra_Map that
+     * sets the parallel partitioning of the vector.
      */
     const Epetra_Map &vector_partitioner () const;
 
     /**
-     *  Output of vector in
-     *  user-defined format in analogy
-     *  to the dealii::Vector<number>
-     *  class.
+     *  Output of vector in user-defined format in analogy to the
+     *  dealii::Vector<number> class.
      */
     void print (const char *format = 0) const;
 
     /**
-     * Print to a stream. @p
-     * precision denotes the desired
-     * precision with which values
-     * shall be printed, @p
-     * scientific whether scientific
-     * notation shall be used. If @p
-     * across is @p true then the
-     * vector is printed in a line,
-     * while if @p false then the
-     * elements are printed on a
+     * Print to a stream. @p precision denotes the desired precision with
+     * which values shall be printed, @p scientific whether scientific
+     * notation shall be used. If @p across is @p true then the vector is
+     * printed in a line, while if @p false then the elements are printed on a
      * separate line each.
      */
     void print (std::ostream       &out,
@@ -1019,42 +805,27 @@ namespace TrilinosWrappers
                 const bool          across     = true) const;
 
     /**
-     * Swap the contents of this
-     * vector and the other vector @p
-     * v. One could do this operation
-     * with a temporary variable and
-     * copying over the data
-     * elements, but this function is
-     * significantly more efficient
-     * since it only swaps the
-     * pointers to the data of the
-     * two vectors and therefore does
-     * not need to allocate temporary
-     * storage and move data
-     * around. Note that the vectors
-     * need to be of the same size
-     * and base on the same map.
+     * Swap the contents of this vector and the other vector @p v. One could
+     * do this operation with a temporary variable and copying over the data
+     * elements, but this function is significantly more efficient since it
+     * only swaps the pointers to the data of the two vectors and therefore
+     * does not need to allocate temporary storage and move data around. Note
+     * that the vectors need to be of the same size and base on the same map.
      *
-     * This function is analog to the
-     * the @p swap function of all C
-     * standard containers. Also,
-     * there is a global function
-     * <tt>swap(u,v)</tt> that simply
-     * calls <tt>u.swap(v)</tt>,
-     * again in analogy to standard
+     * This function is analog to the the @p swap function of all C standard
+     * containers. Also, there is a global function <tt>swap(u,v)</tt> that
+     * simply calls <tt>u.swap(v)</tt>, again in analogy to standard
      * functions.
      */
     void swap (VectorBase &v);
 
     /**
-     * Estimate for the memory
-     * consumption in bytes.
+     * Estimate for the memory consumption in bytes.
      */
     std::size_t memory_consumption () const;
 
     /**
-     * Return a reference to the MPI
-     * communicator object in use with this
+     * Return a reference to the MPI communicator object in use with this
      * object.
      */
     const MPI_Comm &get_mpi_communicator () const;
@@ -1091,57 +862,45 @@ namespace TrilinosWrappers
 
   private:
     /**
-     * Trilinos doesn't allow to
-     * mix additions to matrix
-     * entries and overwriting them
-     * (to make synchronisation of
-     * parallel computations
-     * simpler). The way we do it
-     * is to, for each access
-     * operation, store whether it
-     * is an insertion or an
-     * addition. If the previous
-     * one was of different type,
-     * then we first have to flush
-     * the Trilinos buffers;
-     * otherwise, we can simply go
-     * on.  Luckily, Trilinos has
-     * an object for this which
-     * does already all the
-     * parallel communications in
-     * such a case, so we simply
-     * use their model, which
-     * stores whether the last
-     * operation was an addition or
-     * an insertion.
+     * Trilinos doesn't allow to mix additions to matrix entries and
+     * overwriting them (to make synchronisation of parallel computations
+     * simpler). The way we do it is to, for each access operation, store
+     * whether it is an insertion or an addition. If the previous one was of
+     * different type, then we first have to flush the Trilinos buffers;
+     * otherwise, we can simply go on.  Luckily, Trilinos has an object for
+     * this which does already all the parallel communications in such a case,
+     * so we simply use their model, which stores whether the last operation
+     * was an addition or an insertion.
      */
     Epetra_CombineMode last_action;
 
     /**
-     * A boolean variable to hold
-     * information on whether the
-     * vector is compressed or not.
+     * A boolean variable to hold information on whether the vector is
+     * compressed or not.
      */
     bool compressed;
 
     /**
-     * Whether this vector has ghost elements. This is true
-     * on all processors even if only one of them has any
-     * ghost elements.
+     * Whether this vector has ghost elements. This is true on all processors
+     * even if only one of them has any ghost elements.
      */
     bool has_ghosts;
 
     /**
-     * An Epetra distibuted vector
-     * type. Requires an existing
-     * Epetra_Map for storing data.
+     * An Epetra distibuted vector type. Requires an existing Epetra_Map for
+     * storing data.
      */
     std_cxx1x::shared_ptr<Epetra_FEVector> vector;
 
+    /**
+     * A vector object in Trilinos to be used for collecting the non-local
+     * elements if the vector was constructed with an additional IndexSet
+     * describing ghost elements.
+     */
+    std_cxx1x::shared_ptr<Epetra_MultiVector> nonlocal_vector;
 
     /**
-     * Make the reference class a
-     * friend.
+     * Make the reference class a friend.
      */
     friend class internal::VectorReference;
     friend class Vector;
@@ -1418,6 +1177,9 @@ namespace TrilinosWrappers
     if (fast == false ||
         vector_partitioner().SameAs(v.vector_partitioner())==false)
       vector.reset (new Epetra_FEVector(*v.vector));
+
+    if (v.nonlocal_vector.get() != 0)
+      nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
   }
 
 
@@ -1440,58 +1202,6 @@ namespace TrilinosWrappers
 
 
 
-  inline
-  void
-  VectorBase::compress (::dealii::VectorOperation::values given_last_action)
-  {
-    //Select which mode to send to
-    //Trilinos. Note that we use last_action
-    //if available and ignore what the user
-    //tells us to detect wrongly mixed
-    //operations. Typically given_last_action
-    //is only used on machines that do not
-    //execute an operation (because they have
-    //no own cells for example).
-    Epetra_CombineMode mode = last_action;
-    if (last_action == Zero)
-      {
-        if (given_last_action==::dealii::VectorOperation::add)
-          mode = Add;
-        else if (given_last_action==::dealii::VectorOperation::insert)
-          mode = Insert;
-      }
-
-#ifdef DEBUG
-#  ifdef DEAL_II_WITH_MPI
-    // check that every process has decided
-    // to use the same mode. This will
-    // otherwise result in undefined
-    // behaviour in the call to
-    // GlobalAssemble().
-    double double_mode = mode;
-    Utilities::MPI::MinMaxAvg result
-      = Utilities::MPI::min_max_avg (double_mode,
-                                     dynamic_cast<const Epetra_MpiComm *>
-                                     (&vector_partitioner().Comm())->GetMpiComm());
-    Assert(result.max-result.min<1e-5,
-           ExcMessage ("Not all processors agree whether the last operation on "
-                       "this vector was an addition or a set operation. This will "
-                       "prevent the compress() operation from succeeding."));
-
-#  endif
-#endif
-
-    // Now pass over the information about
-    // what we did last to the vector.
-    const int ierr = vector->GlobalAssemble(mode);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-    last_action = Zero;
-
-    compressed = true;
-  }
-
-
-
   inline
   void
   VectorBase::compress ()
@@ -1515,6 +1225,9 @@ namespace TrilinosWrappers
 
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
+    if (nonlocal_vector.get() != 0)
+      nonlocal_vector->PutScalar(0.);
+
     return *this;
   }
 
@@ -1574,7 +1287,9 @@ namespace TrilinosWrappers
       {
         const size_type row = indices[i];
         const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
-        if (local_row == -1)
+        if (local_row != -1)
+          (*vector)[0][local_row] = values[i];
+        else
           {
             const int ierr = vector->ReplaceGlobalValues (1,
                                                           (const TrilinosWrappers::types::int_type *)(&row),
@@ -1582,8 +1297,11 @@ namespace TrilinosWrappers
             AssertThrow (ierr == 0, ExcTrilinosError(ierr));
             compressed = false;
           }
-        else
-          (*vector)[0][local_row] = values[i];
+        // in set operation, do not use the pre-allocated vector for nonlocal
+        // entries even if it exists. This is to ensure that we really only
+        // set the elements touched by the set() method and not all contained
+        // in the nonlocal entries vector (there is no way to distinguish them
+        // on the receiving processor)
       }
   }
 
@@ -1642,7 +1360,9 @@ namespace TrilinosWrappers
       {
         const size_type row = indices[i];
         const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
-        if (local_row == -1)
+        if (local_row != -1)
+          (*vector)[0][local_row] += values[i];
+        else if (nonlocal_vector.get() == 0)
           {
             const int ierr = vector->SumIntoGlobalValues (1,
                                                           (const TrilinosWrappers::types::int_type *)(&row),
@@ -1651,7 +1371,17 @@ namespace TrilinosWrappers
             compressed = false;
           }
         else
-          (*vector)[0][local_row] += values[i];
+          {
+            // use pre-allocated vector for non-local entries if it exists for
+            // addition operation
+            const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
+            Assert(my_row != -1, 
+                   ExcMessage("Attempted to write into off-processor vector entry "
+                              "that has not be specified as being writable upon "
+                               "initialization"));
+            (*nonlocal_vector)[0][my_row] += values[i];
+            compressed = false;
+          }
       }
   }
 
index ffc29f351e9ff64363accb6858ead20b6dfb5eb0..5df978520e53dc11b9622ff04507d68f51d3e39a 100644 (file)
@@ -1980,17 +1980,19 @@ namespace internal
 
   namespace partition
   {
-    template<typename Worker, bool blocked=false>
+    template<typename Worker>
     class CellWork : public tbb::task
     {
     public:
       CellWork (const Worker &worker_in,
                 const unsigned int partition_in,
-                const internal::MatrixFreeFunctions::TaskInfo &task_info_in)
+                const internal::MatrixFreeFunctions::TaskInfo &task_info_in,
+                const bool is_blocked_in)
         :
         worker (worker_in),
         partition (partition_in),
-        task_info (task_info_in)
+        task_info (task_info_in),
+        is_blocked (is_blocked_in)
       {};
       tbb::task *execute ()
       {
@@ -1998,7 +2000,7 @@ namespace internal
         (task_info.partition_color_blocks_data[partition],
          task_info.partition_color_blocks_data[partition+1]);
         worker(cell_range);
-        if (blocked==true)
+        if (is_blocked==true)
           dummy->spawn (*dummy);
         return NULL;
       }
@@ -2009,92 +2011,82 @@ namespace internal
       const Worker      &worker;
       const unsigned int partition;
       const internal::MatrixFreeFunctions::TaskInfo &task_info;
+      const bool         is_blocked;
     };
 
 
 
-    template<typename Worker, bool blocked=false>
+    template<typename Worker>
     class PartitionWork : public tbb::task
     {
     public:
       PartitionWork (const Worker &function_in,
                      const unsigned int partition_in,
-                     const internal::MatrixFreeFunctions::TaskInfo &task_info_in)
+                     const internal::MatrixFreeFunctions::TaskInfo &task_info_in,
+                     const bool    is_blocked_in = false)
         :
         function (function_in),
         partition (partition_in),
-        task_info (task_info_in)
+        task_info (task_info_in),
+        is_blocked (is_blocked_in)
       {};
       tbb::task *execute ()
       {
-        if (false)
-          {
-            std::pair<unsigned int, unsigned int> cell_range
-            (task_info.partition_color_blocks_data
-             [task_info.partition_color_blocks_row_index[partition]],
-             task_info.partition_color_blocks_data
-             [task_info.partition_color_blocks_row_index[partition+1]]);
-            function(cell_range);
-          }
-        else
+        tbb::empty_task *root = new( tbb::task::allocate_root() )
+          tbb::empty_task;
+        unsigned int evens = task_info.partition_evens[partition];
+        unsigned int odds  = task_info.partition_odds[partition];
+        unsigned int n_blocked_workers =
+          task_info.partition_n_blocked_workers[partition];
+        unsigned int n_workers = task_info.partition_n_workers[partition];
+        std::vector<CellWork<Worker>*> worker(n_workers);
+        std::vector<CellWork<Worker>*> blocked_worker(n_blocked_workers);
+
+        root->set_ref_count(evens+1);
+        for (unsigned int j=0; j<evens; j++)
           {
-            tbb::empty_task *root = new( tbb::task::allocate_root() )
-            tbb::empty_task;
-            unsigned int evens = task_info.partition_evens[partition];
-            unsigned int odds  = task_info.partition_odds[partition];
-            unsigned int n_blocked_workers =
-              task_info.partition_n_blocked_workers[partition];
-            unsigned int n_workers = task_info.partition_n_workers[partition];
-            std::vector<CellWork<Worker,false>*> worker(n_workers);
-            std::vector<CellWork<Worker,true>*> blocked_worker(n_blocked_workers);
-
-            root->set_ref_count(evens+1);
-            for (unsigned int j=0; j<evens; j++)
+            worker[j] = new(root->allocate_child())
+              CellWork<Worker>(function, task_info.
+                               partition_color_blocks_row_index[partition]+2*j,
+                               task_info, false);
+            if (j>0)
               {
-                worker[j] = new(root->allocate_child())
-                CellWork<Worker,false>(function,task_info.
-                                       partition_color_blocks_row_index
-                                       [partition] + 2*j, task_info);
-                if (j>0)
-                  {
-                    worker[j]->set_ref_count(2);
-                    blocked_worker[j-1]->dummy = new(worker[j]->allocate_child())
-                    tbb::empty_task;
-                    worker[j-1]->spawn(*blocked_worker[j-1]);
-                  }
-                else
-                  worker[j]->set_ref_count(1);
-                if (j<evens-1)
+                worker[j]->set_ref_count(2);
+                blocked_worker[j-1]->dummy = new(worker[j]->allocate_child())
+                  tbb::empty_task;
+                worker[j-1]->spawn(*blocked_worker[j-1]);
+              }
+            else
+              worker[j]->set_ref_count(1);
+            if (j<evens-1)
+              {
+                blocked_worker[j] = new(worker[j]->allocate_child())
+                  CellWork<Worker>(function, task_info.
+                                   partition_color_blocks_row_index
+                                   [partition] + 2*j+1, task_info, true);
+              }
+            else
+              {
+                if (odds==evens)
                   {
-                    blocked_worker[j] = new(worker[j]->allocate_child())
-                    CellWork<Worker,true>(function,task_info.
-                                          partition_color_blocks_row_index
-                                          [partition] + 2*j+1, task_info);
+                    worker[evens] = new(worker[j]->allocate_child())
+                      CellWork<Worker>(function, task_info.
+                                       partition_color_blocks_row_index[partition]+2*j+1,
+                                       task_info, false);
+                    worker[j]->spawn(*worker[evens]);
                   }
                 else
                   {
-                    if (odds==evens)
-                      {
-                        worker[evens] = new(worker[j]->allocate_child())
-                        CellWork<Worker,false>(function,
-                                               task_info.
-                                               partition_color_blocks_row_index
-                                               [partition]+2*j+1,task_info);
-                        worker[j]->spawn(*worker[evens]);
-                      }
-                    else
-                      {
-                        tbb::empty_task *child = new(worker[j]->allocate_child())
-                        tbb::empty_task();
-                        worker[j]->spawn(*child);
-                      }
+                    tbb::empty_task *child = new(worker[j]->allocate_child())
+                      tbb::empty_task();
+                    worker[j]->spawn(*child);
                   }
               }
-
-            root->wait_for_all();
-            root->destroy(*root);
           }
-        if (blocked==true)
+
+        root->wait_for_all();
+        root->destroy(*root);
+        if (is_blocked==true)
           dummy->spawn (*dummy);
         return NULL;
       }
@@ -2105,6 +2097,7 @@ namespace internal
       const Worker  &function;
       const unsigned int partition;
       const internal::MatrixFreeFunctions::TaskInfo &task_info;
+      const bool     is_blocked;
     };
 
   } // end of namespace partition
@@ -2150,17 +2143,19 @@ namespace internal
     };
 
 
-    template<typename Worker, bool blocked=false>
+    template<typename Worker>
     class PartitionWork : public tbb::task
     {
     public:
       PartitionWork (const Worker &worker_in,
                      const unsigned int partition_in,
-                     const internal::MatrixFreeFunctions::TaskInfo &task_info_in)
+                     const internal::MatrixFreeFunctions::TaskInfo &task_info_in,
+                     const bool    is_blocked_in)
         :
         worker (worker_in),
         partition (partition_in),
-        task_info (task_info_in)
+        task_info (task_info_in),
+        is_blocked (is_blocked_in)
       {};
       tbb::task *execute ()
       {
@@ -2168,7 +2163,7 @@ namespace internal
                      upper = task_info.partition_color_blocks_data[partition+1];
         parallel_for(tbb::blocked_range<unsigned int>(lower,upper,1),
                      CellWork<Worker> (worker,task_info));
-        if (blocked==true)
+        if (is_blocked==true)
           dummy->spawn (*dummy);
         return NULL;
       }
@@ -2179,6 +2174,7 @@ namespace internal
       const Worker &worker;
       const unsigned int partition;
       const internal::MatrixFreeFunctions::TaskInfo &task_info;
+      const bool is_blocked;
     };
 
   } // end of namespace color
@@ -2273,9 +2269,9 @@ MatrixFree<dim, Number>::cell_loop
           root->set_ref_count(evens+1);
           unsigned int n_blocked_workers = task_info.n_blocked_workers;
           unsigned int n_workers = task_info.n_workers;
-          std::vector<internal::partition::PartitionWork<Worker,false>*>
+          std::vector<internal::partition::PartitionWork<Worker>*>
           worker(n_workers);
-          std::vector<internal::partition::PartitionWork<Worker,true>*>
+          std::vector<internal::partition::PartitionWork<Worker>*>
           blocked_worker(n_blocked_workers);
           internal::MPIComCompress<OutVector> *worker_compr =
             new(root->allocate_child())
@@ -2286,8 +2282,8 @@ MatrixFree<dim, Number>::cell_loop
               if (j>0)
                 {
                   worker[j] = new(root->allocate_child())
-                  internal::partition::PartitionWork<Worker,false>
-                  (func,2*j,task_info);
+                  internal::partition::PartitionWork<Worker>
+                  (func,2*j,task_info,false);
                   worker[j]->set_ref_count(2);
                   blocked_worker[j-1]->dummy = new(worker[j]->allocate_child())
                   tbb::empty_task;
@@ -2299,8 +2295,8 @@ MatrixFree<dim, Number>::cell_loop
               else
                 {
                   worker[j] = new(worker_compr->allocate_child())
-                  internal::partition::PartitionWork<Worker,false>
-                  (func,2*j,task_info);
+                  internal::partition::PartitionWork<Worker>
+                  (func,2*j,task_info,false);
                   worker[j]->set_ref_count(2);
                   internal::MPIComDistribute<InVector> *worker_dist =
                     new (worker[j]->allocate_child())
@@ -2310,16 +2306,16 @@ MatrixFree<dim, Number>::cell_loop
               if (j<evens-1)
                 {
                   blocked_worker[j] = new(worker[j]->allocate_child())
-                  internal::partition::PartitionWork<Worker,true>
-                  (func,2*j+1,task_info);
+                  internal::partition::PartitionWork<Worker>
+                  (func,2*j+1,task_info,true);
                 }
               else
                 {
                   if (odds==evens)
                     {
                       worker[evens] = new(worker[j]->allocate_child())
-                      internal::partition::PartitionWork<Worker,false>
-                      (func,2*j+1,task_info);
+                      internal::partition::PartitionWork<Worker>
+                      (func,2*j+1,task_info,false);
                       worker[j]->spawn(*worker[evens]);
                     }
                   else
@@ -2348,8 +2344,8 @@ MatrixFree<dim, Number>::cell_loop
               unsigned int n_blocked_workers = odds-(odds+evens+1)%2;
               unsigned int n_workers = task_info.partition_color_blocks_data.size()-1-
                                        n_blocked_workers;
-              std::vector<internal::color::PartitionWork<Worker,false>*> worker(n_workers);
-              std::vector<internal::color::PartitionWork<Worker,true>*> blocked_worker(n_blocked_workers);
+              std::vector<internal::color::PartitionWork<Worker>*> worker(n_workers);
+              std::vector<internal::color::PartitionWork<Worker>*> blocked_worker(n_blocked_workers);
               unsigned int worker_index = 0, slice_index = 0;
               unsigned int spawn_index =  0, spawn_index_new = 0;
               int spawn_index_child = -2;
@@ -2362,10 +2358,10 @@ MatrixFree<dim, Number>::cell_loop
                   spawn_index_new = worker_index;
                   if (part == 0)
                     worker[worker_index] = new(worker_compr->allocate_child())
-                    internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
+                    internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                   else
                     worker[worker_index] = new(root->allocate_child())
-                    internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
+                    internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                   slice_index++;
                   for (; slice_index<task_info.partition_color_blocks_row_index[part+1];
                        slice_index++)
@@ -2373,7 +2369,7 @@ MatrixFree<dim, Number>::cell_loop
                       worker[worker_index]->set_ref_count(1);
                       worker_index++;
                       worker[worker_index] = new (worker[worker_index-1]->allocate_child())
-                      internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
+                      internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                     }
                   worker[worker_index]->set_ref_count(2);
                   if (part>0)
@@ -2402,14 +2398,14 @@ MatrixFree<dim, Number>::cell_loop
                       if (part<task_info.partition_color_blocks_row_index.size()-2)
                         {
                           blocked_worker[part/2] = new(worker[worker_index-1]->allocate_child())
-                          internal::color::PartitionWork<Worker,true>(func,slice_index,task_info);
+                          internal::color::PartitionWork<Worker>(func,slice_index,task_info,true);
                           slice_index++;
                           if (slice_index<
                               task_info.partition_color_blocks_row_index[part+1])
                             {
                               blocked_worker[part/2]->set_ref_count(1);
                               worker[worker_index] = new(blocked_worker[part/2]->allocate_child())
-                              internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
+                              internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                               slice_index++;
                             }
                           else
@@ -2428,7 +2424,7 @@ MatrixFree<dim, Number>::cell_loop
                               worker_index++;
                             }
                           worker[worker_index] = new (worker[worker_index-1]->allocate_child())
-                          internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
+                          internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                         }
                       spawn_index_child = worker_index;
                       worker_index++;
index 402f7822ea1751240586290ddb213f47c8de1c1d..7d0f7f62f8f4350acb3ffa83364d21f0460461f2 100644 (file)
@@ -172,6 +172,10 @@ namespace MeshWorker
     /// The block structure of the system
     SmartPointer<const BlockInfo,DoFInfo<dim,spacedim> > block_info;
 
+    /**
+     * The structure refers to a cell with level data instead of
+     * active data.
+     */
     bool level_cell;
   private:
     /**
@@ -280,6 +284,12 @@ namespace MeshWorker
      * face is available.
      */
     bool exterior_face_available[GeometryInfo<dim>::faces_per_cell];
+
+    /**
+     * A flag to specify if the current object has been set to a valid
+     * cell.
+     */
+    bool cell_valid;
   };
 
 //----------------------------------------------------------------------//
@@ -417,7 +427,7 @@ namespace MeshWorker
   inline
   DoFInfoBox<dim, DOFINFO>::DoFInfoBox(const DOFINFO &seed)
     :
-    cell(seed)
+    cell(seed), cell_valid(true)
   {
     for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
       {
@@ -433,7 +443,7 @@ namespace MeshWorker
   inline
   DoFInfoBox<dim, DOFINFO>::DoFInfoBox(const DoFInfoBox<dim, DOFINFO> &other)
     :
-    cell(other.cell)
+    cell(other.cell), cell_valid(other.cell_valid)
   {
     for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
       {
@@ -449,6 +459,7 @@ namespace MeshWorker
   inline void
   DoFInfoBox<dim, DOFINFO>::reset ()
   {
+      cell_valid = false;
     for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
       {
         interior_face_available[i] = false;
@@ -462,6 +473,9 @@ namespace MeshWorker
   inline void
   DoFInfoBox<dim, DOFINFO>::assemble (ASSEMBLER &assembler) const
   {
+    if (!cell_valid)
+      return;
+
     assembler.assemble(cell);
     for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
       {
index 1bd8c63b5c4642674fcaf5dc96b953f91df385b7..4eb299711e6c3e58c98b0a470eb2dd920cdcde9d 100644 (file)
@@ -386,7 +386,7 @@ namespace MeshWorker
      * these grid entities should be performed.
      *
      * If the parameter <tt>force</tt> is true, then all quadrature
-     * sets are filled with new quadrature ruels. If it is false, then
+     * sets are filled with new quadrature rules. If it is false, then
      * only empty rules are changed.
      */
     void initialize_gauss_quadrature(unsigned int n_cell_points,
index 6595f91d94dc570a8b3b46482ec7b33a8b916520..22110233a2c0d1c9891742a42f1d26ca7c31bc8c 100644 (file)
@@ -81,7 +81,8 @@ namespace MeshWorker
        */
       LoopControl()
       : own_cells(true), ghost_cells(false),
-        faces_to_ghost(LoopControl::one), own_faces(LoopControl::one)
+        faces_to_ghost(LoopControl::one), own_faces(LoopControl::one),
+    cells_first(true)
       {
       }
 
@@ -117,80 +118,19 @@ namespace MeshWorker
        * Loop over faces between two locally owned cells:
        * - never: do not assemble face terms
        * - one: assemble once (always coming from the finer side)
-       * - both: assemble each face twice           ALSO IF FINER?!?!?
-       * Default is one_side.
+       * - both: assemble each face twice (not implemented for hanging nodes!)
+       * Default is one.
        */
       FaceOption own_faces;
 
 
-
-      /**
-       * Based on the flags in this class, decide if this face needs to be
-       * assembled.
-       */
-      template <class ITERATOR>
-      bool operator() (const ITERATOR& cell, unsigned int face, bool is_level) const;
       /**
-       * Based on the flags in this class, decide if this cell needs to be
-       * assembled.
+       * Flag to determine if cells integrals should be done before or after
+       * face integrals. Default is t
        */
-      template <class ITERATOR>
-      bool operator() (const ITERATOR& cell, bool is_level) const;
+      bool cells_first;
   };
   
-  
-  template <class ITERATOR>
-  bool LoopControl::operator() (const ITERATOR& cell, unsigned int face, bool is_level) const
-  {
-    const ITERATOR neighbor = cell->neighbor(face);
-    const bool c_local = (is_level)
-                        ? (cell->is_locally_owned())
-                        : (cell->is_locally_owned_on_level());
-    const bool n_local = (is_level)
-                        ? (neighbor->is_locally_owned())
-                        : (neighbor->is_locally_owned_on_level());
-    
-    if (!c_local && !n_local)
-      return false;
-
-    if (c_local && n_local)
-      {
-        if (own_faces==LoopControl::never)
-          return false;
-
-        //TODO:
-//        if (cell->neighbor_is_coarser(face_no))
-//          return false;
-
-        if (own_faces==LoopControl::one && neighbor < cell)
-          return false;
-
-        Assert(own_faces==LoopControl::both, ExcInternalError());
-        return true;
-      }
-    else
-      {
-        // interface between owned and ghost cell
-
-        // TODO
-
-
-      }
-  }
-  template <class ITERATOR>
-  bool LoopControl::operator() (const ITERATOR &cell, bool is_level) const
-  {
-    const bool c_local = (is_level)
-                        ? (cell->is_locally_owned())
-                        : (cell->is_locally_owned_on_level());
-    if (own_cells && c_local)
-      return true;
-    if (ghost_cells && cell->is_ghost())
-      return true;
-    return false;
-  }
-
-
 
 
   /**
@@ -220,10 +160,7 @@ namespace MeshWorker
    *        integrals are to be  dealt with first. Note that independent of the
    *        value of this flag, cell and face integrals of a given cell are
    *        all taken care of before moving to the next cell.
-   * @param unique_faces_only determines, that a face between two cells
-   * of the same level is processed only from the cell which is less
-   * than its neighbor. If this parameter is <tt>false</tt> these faces
-   * are processed from both cells.
+   * @param loop_control control structure to specify what actions should be performed.
    *
    * @ingroup MeshWorker
    * @author Guido Kanschat
@@ -239,7 +176,6 @@ namespace MeshWorker
     const std_cxx1x::function<void (DOFINFO &, DOFINFO &,
                                     typename INFOBOX::CellInfo &,
                                     typename INFOBOX::CellInfo &)> &face_worker,
-    const bool cells_first,
     const LoopControl & loop_control)
   {
     const bool ignore_subdomain = (cell->get_triangulation().locally_owned_subdomain()
@@ -249,25 +185,27 @@ namespace MeshWorker
                                ? cell->level_subdomain_id()
                                : cell->subdomain_id();
      
+     const bool own_cell = ignore_subdomain || (csid == cell->get_triangulation().locally_owned_subdomain());
+
+     dof_info.reset();
+
      if ((!ignore_subdomain) && (csid == numbers::artificial_subdomain_id))
        return;
-     
-     dof_info.reset();
+
      dof_info.cell.reinit(cell);
+     dof_info.cell_valid = true;
      
     const bool integrate_cell          = (cell_worker != 0);
     const bool integrate_boundary      = (boundary_worker != 0);
     const bool integrate_interior_face = (face_worker != 0);
 
-    dof_info.reset();
-
-    dof_info.cell.reinit(cell);
     if (integrate_cell)
       info.cell.reinit(dof_info.cell);
     // Execute this, if cells
     // have to be dealt with
     // before faces
-    if (integrate_cell && cells_first)
+    if (integrate_cell && loop_control.cells_first &&
+        ((loop_control.own_cells && own_cell) || (loop_control.ghost_cells && !own_cell)))
       cell_worker(dof_info.cell, info.cell);
 
     // Call the callback function in
@@ -282,7 +220,8 @@ namespace MeshWorker
           typename ITERATOR::AccessorType::Container::face_iterator face = cell->face(face_no);
           if (cell->at_boundary(face_no))
             {
-              if (integrate_boundary)
+              // only integrate boundary faces of own cells
+              if (integrate_boundary && own_cell)
                 {
                   dof_info.interior_face_available[face_no] = true;
                   dof_info.interior[face_no].reinit(cell, face, face_no);
@@ -295,6 +234,28 @@ namespace MeshWorker
               // Interior face
               TriaIterator<typename ITERATOR::AccessorType> neighbor = cell->neighbor(face_no);
 
+              types::subdomain_id neighbid = numbers::artificial_subdomain_id;
+              if (neighbor->is_level_cell())
+                neighbid = neighbor->level_subdomain_id();
+              //subdomain id is only valid for active cells
+              else if (neighbor->active())
+                neighbid = neighbor->subdomain_id();
+
+              const bool own_neighbor = ignore_subdomain ||
+                  (neighbid == cell->get_triangulation().locally_owned_subdomain());
+
+              // skip all faces between two ghost cells
+              if (!own_cell && !own_neighbor)
+                continue;
+
+              // skip if the user doesn't want faces between own cells
+              if (own_cell && own_neighbor && loop_control.own_faces==LoopControl::never)
+                continue;
+
+              // skip face to ghost
+              if (own_cell != own_neighbor && loop_control.faces_to_ghost==LoopControl::never)
+                              continue;
+
               // Deal with
               // refinement edges
               // from the refined
@@ -310,6 +271,12 @@ namespace MeshWorker
                   Assert(!cell->has_children(), ExcInternalError());
                   Assert(!neighbor->has_children(), ExcInternalError());
 
+                  // skip if only one processor needs to assemble the face
+                  // to a ghost cell and the fine cell is not ours.
+                  if (!own_cell
+                      && loop_control.faces_to_ghost == LoopControl::one)
+                    continue;
+
                   const std::pair<unsigned int, unsigned int> neighbor_face_no
                     = cell->neighbor_of_coarser_neighbor(face_no);
                   const typename ITERATOR::AccessorType::Container::face_iterator nface
@@ -328,20 +295,38 @@ namespace MeshWorker
                 }
               else
                 {
-                  // Neighbor is
-                  // on same
-                  // level, but
-                  // only do this
-                  // from one side.
-                  if (loop_control.own_faces != LoopControl::both && (neighbor < cell)) continue;
-
-                  // If iterator
-                  // is active
-                  // and neighbor
-                  // is refined,
-                  // skip
+                  // If iterator is active and neighbor is refined, skip
                   // internal face.
                   if (internal::is_active_iterator(cell) && neighbor->has_children())
+                    {
+                      Assert(loop_control.own_faces != LoopControl::both, ExcMessage(
+                          "Assembling from both sides for own_faces is not "
+                          "supported with hanging nodes!"));
+                      continue;
+                    }
+
+                  // Now neighbor is on same level, double-check this:
+                  Assert(cell->level()==neighbor->level(), ExcInternalError());
+
+                  // only do faces on same level from one side (unless
+                  // LoopControl says otherwise)
+                  if (own_cell && own_neighbor
+                      && loop_control.own_faces == LoopControl::one
+                      && (neighbor < cell))
+                    continue;
+
+                  // independent of loop_control.faces_to_ghost,
+                  // we only look at faces to ghost on the same level once
+                  // (only where own_cell=true and own_neighbor=false)
+                  if (!own_cell)
+                    continue;
+
+                  // now only one processor assembles faces_to_ghost. This
+                  // logic is based on the subdomain id and is handled inside
+                  // operator<.
+                  if (own_cell && !own_neighbor
+                      && loop_control.faces_to_ghost == LoopControl::one
+                      && (neighbor < cell))
                     continue;
 
                   const unsigned int neighbor_face_no = cell->neighbor_face_no(face_no);
@@ -368,7 +353,8 @@ namespace MeshWorker
 
     // Execute this, if faces
     // have to be handled first
-    if (integrate_cell && !cells_first)
+    if (integrate_cell && !loop_control.cells_first &&
+      ((loop_control.own_cells && own_cell) || (loop_control.ghost_cells && !own_cell)))
       cell_worker(dof_info.cell, info.cell);
   }
 
@@ -389,6 +375,48 @@ namespace MeshWorker
    * @ingroup MeshWorker
    * @author Guido Kanschat, 2009
    */
+  template<int dim, int spacedim, class DOFINFO, class INFOBOX, class ASSEMBLER, class ITERATOR>
+    void loop(ITERATOR begin,
+              typename identity<ITERATOR>::type end,
+              DOFINFO &dinfo,
+              INFOBOX &info,
+              const std_cxx1x::function<void (DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker,
+              const std_cxx1x::function<void (DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker,
+              const std_cxx1x::function<void (DOFINFO &, DOFINFO &,
+                                              typename INFOBOX::CellInfo &,
+                                              typename INFOBOX::CellInfo &)> &face_worker,
+              ASSEMBLER &assembler,
+              const LoopControl &lctrl = LoopControl())
+    {
+      DoFInfoBox<dim, DOFINFO> dof_info(dinfo);
+
+      assembler.initialize_info(dof_info.cell, false);
+      for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+        {
+          assembler.initialize_info(dof_info.interior[i], true);
+          assembler.initialize_info(dof_info.exterior[i], true);
+        }
+
+      // Loop over all cells
+  #ifdef DEAL_II_MESHWORKER_PARALLEL
+      WorkStream::run(begin, end,
+                      std_cxx1x::bind(&cell_action<INFOBOX, DOFINFO, dim, spacedim, ITERATOR>,
+                                      std_cxx1x::_1, std_cxx1x::_3, std_cxx1x::_2,
+                                      cell_worker, boundary_worker, face_worker, lctrl),
+                      std_cxx1x::bind(&internal::assemble<dim,DOFINFO,ASSEMBLER>, std_cxx1x::_1, &assembler),
+                      info, dof_info);
+  #else
+      for (ITERATOR cell = begin; cell != end; ++cell)
+        {
+          cell_action<INFOBOX,DOFINFO,dim,spacedim>(cell, dof_info,
+                                                    info, cell_worker,
+                                                    boundary_worker, face_worker,
+                                                    lctrl);
+          dof_info.assemble(assembler);
+        }
+  #endif
+    }
+
   template<int dim, int spacedim, class DOFINFO, class INFOBOX, class ASSEMBLER, class ITERATOR>
   void loop(ITERATOR begin,
             typename identity<ITERATOR>::type end,
@@ -400,40 +428,28 @@ namespace MeshWorker
                                             typename INFOBOX::CellInfo &,
                                             typename INFOBOX::CellInfo &)> &face_worker,
             ASSEMBLER &assembler,
-            bool cells_first = true,
-            bool unique_faces_only = true)
-  {
-    DoFInfoBox<dim, DOFINFO> dof_info(dinfo);
+            bool cells_first,
+            bool unique_faces_only = true) DEAL_II_DEPRECATED;
 
-    assembler.initialize_info(dof_info.cell, false);
-    for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
-      {
-        assembler.initialize_info(dof_info.interior[i], true);
-        assembler.initialize_info(dof_info.exterior[i], true);
-      }
+  template<int dim, int spacedim, class DOFINFO, class INFOBOX, class ASSEMBLER, class ITERATOR>
+  void loop(ITERATOR begin,
+              typename identity<ITERATOR>::type end,
+              DOFINFO &dinfo,
+              INFOBOX &info,
+              const std_cxx1x::function<void (DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker,
+              const std_cxx1x::function<void (DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker,
+              const std_cxx1x::function<void (DOFINFO &, DOFINFO &,
+                                              typename INFOBOX::CellInfo &,
+                                              typename INFOBOX::CellInfo &)> &face_worker,
+              ASSEMBLER &assembler,
+              bool cells_first,
+              bool unique_faces_only)
+  {
+      LoopControl lctrl;
+      lctrl.cells_first = cells_first;
+      lctrl.own_faces = (unique_faces_only)?LoopControl::one:LoopControl::both;
 
-    LoopControl lctrl;
-    lctrl.own_faces = (unique_faces_only)?LoopControl::one:LoopControl::both;
-
-    // Loop over all cells
-#ifdef DEAL_II_MESHWORKER_PARALLEL
-    WorkStream::run(begin, end,
-                    std_cxx1x::bind(&cell_action<INFOBOX, DOFINFO, dim, spacedim, ITERATOR>,
-                                    std_cxx1x::_1, std_cxx1x::_3, std_cxx1x::_2,
-                                    cell_worker, boundary_worker, face_worker, cells_first, lctrl),
-                    std_cxx1x::bind(&internal::assemble<dim,DOFINFO,ASSEMBLER>, std_cxx1x::_1, &assembler),
-                    info, dof_info);
-#else
-    for (ITERATOR cell = begin; cell != end; ++cell)
-      {
-        cell_action<INFOBOX,DOFINFO,dim,spacedim>(cell, dof_info,
-                                                  info, cell_worker,
-                                                  boundary_worker, face_worker,
-                                                  cells_first,
-                                                  lctrl);
-        dof_info.assemble(assembler);
-      }
-#endif
+      loop<dim,spacedim>(begin, end, dinfo, info, cell_worker, boundary_worker, face_worker, assembler, lctrl);
   }
 
   /**
@@ -456,7 +472,7 @@ namespace MeshWorker
                                                         IntegrationInfo<dim, spacedim> &,
                                                         IntegrationInfo<dim, spacedim> &)> &face_worker,
                         ASSEMBLER &assembler,
-                        bool cells_first = true) DEAL_II_DEPRECATED;
+                        bool cells_first) DEAL_II_DEPRECATED;
 
 
   template<int dim, int spacedim, class ITERATOR, class ASSEMBLER>
@@ -501,7 +517,17 @@ namespace MeshWorker
                         IntegrationInfoBox<dim, spacedim> &box,
                         const LocalIntegrator<dim, spacedim> &integrator,
                         ASSEMBLER &assembler,
-                        bool cells_first = true)
+                        bool cells_first)
+DEAL_II_DEPRECATED;
+
+  template<int dim, int spacedim, class ITERATOR, class ASSEMBLER>
+  void integration_loop(ITERATOR begin,
+                        typename identity<ITERATOR>::type end,
+                        DoFInfo<dim, spacedim> &dof_info,
+                        IntegrationInfoBox<dim, spacedim> &box,
+                        const LocalIntegrator<dim, spacedim> &integrator,
+                        ASSEMBLER &assembler,
+                        bool cells_first)
   {
     std_cxx1x::function<void (DoFInfo<dim>&, IntegrationInfo<dim, spacedim>&)> cell_worker;
     std_cxx1x::function<void (DoFInfo<dim>&, IntegrationInfo<dim, spacedim>&)> boundary_worker;
@@ -526,7 +552,40 @@ namespace MeshWorker
      cells_first);
   }
 
+  /**
+   * As above but using LoopControl
+   */
+  template<int dim, int spacedim, class ITERATOR, class ASSEMBLER>
+  void integration_loop(ITERATOR begin,
+                        typename identity<ITERATOR>::type end,
+                        DoFInfo<dim, spacedim> &dof_info,
+                        IntegrationInfoBox<dim, spacedim> &box,
+                        const LocalIntegrator<dim, spacedim> &integrator,
+                        ASSEMBLER &assembler,
+                        const LoopControl &lctrl = LoopControl())
+  {
+    std_cxx1x::function<void (DoFInfo<dim>&, IntegrationInfo<dim, spacedim>&)> cell_worker;
+    std_cxx1x::function<void (DoFInfo<dim>&, IntegrationInfo<dim, spacedim>&)> boundary_worker;
+    std_cxx1x::function<void (DoFInfo<dim> &, DoFInfo<dim> &,
+                              IntegrationInfo<dim, spacedim> &,
+                              IntegrationInfo<dim, spacedim> &)> face_worker;
+    if (integrator.use_cell)
+      cell_worker = std_cxx1x::bind(&LocalIntegrator<dim, spacedim>::cell, &integrator, std_cxx1x::_1, std_cxx1x::_2);
+    if (integrator.use_boundary)
+      boundary_worker = std_cxx1x::bind(&LocalIntegrator<dim, spacedim>::boundary, &integrator, std_cxx1x::_1, std_cxx1x::_2);
+    if (integrator.use_face)
+      face_worker = std_cxx1x::bind(&LocalIntegrator<dim, spacedim>::face, &integrator, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3, std_cxx1x::_4);
 
+    loop<dim, spacedim>
+    (begin, end,
+     dof_info,
+     box,
+     cell_worker,
+     boundary_worker,
+     face_worker,
+     assembler,
+     lctrl);
+  }
 
 }
 
index 9f279b3c8bc7f41b81d32dbf0bfbc9f14e850315..224729da485d29a235ddda1cf043a6efe7b5c08c 100644 (file)
@@ -824,24 +824,13 @@ namespace MeshWorker
     {
       AssertDimension(M.m(), i1.size());
       AssertDimension(M.n(), i2.size());
+      Assert(mg_constrained_dofs == 0, ExcInternalError());
+//TODO: Possibly remove this function all together
 
-      if (mg_constrained_dofs == 0)
-        {
-          for (unsigned int j=0; j<i1.size(); ++j)
-            for (unsigned int k=0; k<i2.size(); ++k)
-              if (std::fabs(M(j,k)) >= threshold)
-                G.add(i1[j], i2[k], M(j,k));
-        }
-      else
-        {
-          for (unsigned int j=0; j<i1.size(); ++j)
-            for (unsigned int k=0; k<i2.size(); ++k)
-              if (std::fabs(M(j,k)) >= threshold)
-                {
-                  if (!mg_constrained_dofs->continuity_across_refinement_edges())
-                    G.add(i1[j], i2[k], M(j,k));
-                }
-        }
+      for (unsigned int j=0; j<i1.size(); ++j)
+        for (unsigned int k=0; k<i2.size(); ++k)
+          if (std::fabs(M(j,k)) >= threshold)
+            G.add(i1[j], i2[k], M(j,k));
     }
 
 
@@ -868,28 +857,35 @@ namespace MeshWorker
         {
           for (unsigned int j=0; j<i1.size(); ++j)
             for (unsigned int k=0; k<i2.size(); ++k)
-              if (std::fabs(M(j,k)) >= threshold)
-                if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
-                    !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
-                  {
-                    if (mg_constrained_dofs->set_boundary_values())
-                      {
-                        // At the boundary, only enter the term on the
-                        // diagonal, but not the coupling terms
-                        if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
-                             !mg_constrained_dofs->is_boundary_index(level, i2[k]))
-                            ||
-                            (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
-                             mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
-                             i1[j] == i2[k]))
-                          G.add(i1[j], i2[k], M(j,k));
-                      }
-                    else
-                      G.add(i1[j], i2[k], M(j,k));
-                  }
-        }
+             {
+               // Only enter the local values into the global matrix,
+               //  if the value is larger than the threshold
+               if (std::fabs(M(j,k)) < threshold)
+                 continue;
+               
+               // Do not enter, if either the row or the column
+               // corresponds to an index on the refinement edge. The
+               // level problems are solved with homogeneous
+               // Dirichlet boundary conditions, therefore we
+               // eliminate these rows and columns. The corresponding
+               // matrix entries are entered by assemble_in() and
+               // assemble_out().
+                if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) ||
+                    mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+                 continue;
+               
+               // At the boundary, only enter the term on the
+               // diagonal, but not the coupling terms
+               if ((mg_constrained_dofs->is_boundary_index(level, i1[j]) ||
+                    mg_constrained_dofs->is_boundary_index(level, i2[k])) &&
+                   (i1[j] != i2[k]))
+                 continue;
+               
+               G.add(i1[j], i2[k], M(j,k));
+             }
+       }
     }
-
+    
 
     template <class MATRIX>
     inline void
@@ -1085,28 +1081,18 @@ namespace MeshWorker
         {
           if (level1 == level2)
             {
-              if (mg_constrained_dofs == 0)
-                {
-                  assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
-                  assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices);
-                  assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices);
-                  assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
-                }
-              else
-                {
-                  assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
-                  assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
-                  assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
-                  assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
-                }
-            }
+             assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
+             assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
+             assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
+             assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+           }
           else
             {
               Assert(level1 > level2, ExcInternalError());
               // Do not add info2.M1,
               // which is done by
               // the coarser cell
-              assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
+              assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
               if (level1>0)
                 {
                   assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1);
@@ -1122,20 +1108,10 @@ namespace MeshWorker
 
             if (level1 == level2)
               {
-                if (mg_constrained_dofs == 0)
-                  {
-                    assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
-                    assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column]);
-                    assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column]);
-                    assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column]);
-                  }
-                else
-                  {
-                    assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
-                    assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1);
-                    assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1);
-                    assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
-                  }
+               assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
+               assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1);
+               assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1);
+               assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
               }
             else
               {
@@ -1143,7 +1119,7 @@ namespace MeshWorker
                 // Do not add info2.M1,
                 // which is done by
                 // the coarser cell
-                assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
+                assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
                 if (level1>0)
                   {
                     assemble_up((*flux_up)[level1],info1.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
index 610e0344a49bcaf2b909f9aef3da305450949fc8..9d528f1204bf147254ab9b505e614b36bb503ac9 100644 (file)
@@ -106,47 +106,53 @@ public:
    * Return the indices of dofs for each level that lie on the
    * boundary of the domain.
    */
-// TODO: remove
   const std::vector<std::set<types::global_dof_index> > &
   get_boundary_indices () const;
 
   /**
-   * Return the indices of dofs for each level that lie on the
-   * boundary of the domain.
-   */
-  // TODO: remove
-  const std::vector<std::set<types::global_dof_index> > &
-  get_non_refinement_edge_indices () const;
-
-  /**
+   * @deprecated Use at_refinement_edge() if possible, else
+   * get_refinement_edge_indices(unsigned int).
+   *
    * Return the indices of dofs for each level that lie on the
    * refinement edge (i.e. are on faces between cells of this level
    * and cells on the level below).
    */
-  // TODO: remove
   const std::vector<std::vector<bool> > &
-  get_refinement_edge_indices () const;
+  get_refinement_edge_indices () const DEAL_II_DEPRECATED;
 
   /**
+   * Return the indices of dofs on the given level that lie on an
+   * refinement edge (dofs on faces to neighbors that are coarser)
+   */
+  const IndexSet &
+  get_refinement_edge_indices (unsigned int level) const;
+
+  /**
+   * @deprecated Use at_refinement_edge_boundary() if possible, else
+   * use get_refinement_edge_boundary_indices().
+   *
    * Return the indices of dofs for each level that are in the
    * intersection of the sets returned by get_boundary_indices() and
    * get_refinement_edge_indices().
    */
-  // TODO: remove
   const std::vector<std::vector<bool> > &
-  get_refinement_edge_boundary_indices () const;
+  get_refinement_edge_boundary_indices () const DEAL_II_DEPRECATED;
 
   /**
-   * Return if boundary_indices need to be set or not.
+   * Return indices of all dofs that are on boundary faces on the given level 
+   * if the cell has refinement edge indices (i.e. has a coarser neighbor).
    */
-
-  bool set_boundary_values () const;
+  const IndexSet &
+  get_refinement_edge_boundary_indices (unsigned int level) const;
 
   /**
-   * Return if the finite element requires continuity across
-   * refinement edges.
+   * @deprecated The function is_boundary_index() now returns false if
+   * no boundary values are set.
+   *
+   * Return if boundary_indices need to be set or not.
    */
-  bool continuity_across_refinement_edges () const;
+  bool set_boundary_values () const DEAL_II_DEPRECATED;
+
 private:
 
   /**
@@ -154,17 +160,11 @@ private:
    */
   std::vector<std::set<types::global_dof_index> > boundary_indices;
 
-  /**
-   * The degrees of freedom on egdges that are not a refinement edge
-   * between a level and coarser cells.
-   */
-  std::vector<std::set<types::global_dof_index> > non_refinement_edge_indices;
-
   /**
    * The degrees of freedom on the refinement edge between a level and
    * coarser cells.
    */
-  std::vector<std::vector<bool> > refinement_edge_indices;
+  std::vector<IndexSet> refinement_edge_indices;
 
   /**
    * The degrees of freedom on the refinement edge between a level and
@@ -172,7 +172,17 @@ private:
    *
    * This is a subset of #refinement_edge_indices.
    */
-  std::vector<std::vector<bool> > refinement_edge_boundary_indices;
+  std::vector<IndexSet> refinement_edge_boundary_indices;
+
+  /**
+   * old data structure only filled on demand
+   */
+  mutable std::vector<std::vector<bool> > refinement_edge_boundary_indices_old;
+
+  /**
+   * old data structure only filled on demand
+   */
+  mutable std::vector<std::vector<bool> > refinement_edge_indices_old;
 };
 
 
@@ -184,16 +194,15 @@ MGConstrainedDoFs::initialize(const DoFHandler<dim,spacedim> &dof)
   const unsigned int nlevels = dof.get_tria().n_global_levels();
   refinement_edge_indices.resize(nlevels);
   refinement_edge_boundary_indices.resize(nlevels);
-  non_refinement_edge_indices.resize(nlevels);
+  refinement_edge_indices_old.clear();
+  refinement_edge_boundary_indices_old.clear();
   for (unsigned int l=0; l<nlevels; ++l)
     {
-      refinement_edge_indices[l].resize(dof.n_dofs(l));
-      refinement_edge_boundary_indices[l].resize(dof.n_dofs(l));
+      refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
+      refinement_edge_boundary_indices[l] = IndexSet(dof.n_dofs(l));
     }
   MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices,
                                          refinement_edge_boundary_indices);
-
-  MGTools::extract_non_interface_dofs (dof, non_refinement_edge_indices);
 }
 
 
@@ -209,19 +218,19 @@ MGConstrainedDoFs::initialize(
   boundary_indices.resize(nlevels);
   refinement_edge_indices.resize(nlevels);
   refinement_edge_boundary_indices.resize(nlevels);
-  non_refinement_edge_indices.resize(nlevels);
+  refinement_edge_indices_old.clear();
+  refinement_edge_boundary_indices_old.clear();
 
   for (unsigned int l=0; l<nlevels; ++l)
     {
       boundary_indices[l].clear();
-      refinement_edge_indices[l].resize(dof.n_dofs(l));
-      refinement_edge_boundary_indices[l].resize(dof.n_dofs(l));
+      refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
+      refinement_edge_boundary_indices[l] = IndexSet(dof.n_dofs(l));
     }
 
   MGTools::make_boundary_list (dof, function_map, boundary_indices, component_mask);
   MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices,
                                          refinement_edge_boundary_indices);
-  MGTools::extract_non_interface_dofs (dof, non_refinement_edge_indices);
 }
 
 
@@ -238,8 +247,8 @@ MGConstrainedDoFs::clear()
   for (unsigned int l=0; l<refinement_edge_boundary_indices.size(); ++l)
     refinement_edge_boundary_indices[l].clear();
 
-  for (unsigned int l=0; l<non_refinement_edge_indices.size(); ++l)
-    non_refinement_edge_indices[l].clear();
+  refinement_edge_indices_old.clear();
+  refinement_edge_boundary_indices_old.clear();
 }
 
 
@@ -248,11 +257,11 @@ bool
 MGConstrainedDoFs::is_boundary_index (const unsigned int level,
                                       const types::global_dof_index index) const
 {
-  AssertIndexRange(level, boundary_indices.size());
-  if (boundary_indices[level].find(index) != boundary_indices[level].end())
-    return true;
-  else
+  if (boundary_indices.size() == 0)
     return false;
+  
+  AssertIndexRange(level, boundary_indices.size());
+  return (boundary_indices[level].find(index) != boundary_indices[level].end());
 }
 
 inline
@@ -260,12 +269,7 @@ bool
 MGConstrainedDoFs::non_refinement_edge_index (const unsigned int level,
                                               const types::global_dof_index index) const
 {
-  AssertIndexRange(level, non_refinement_edge_indices.size());
-
-  if (non_refinement_edge_indices[level].find(index) != non_refinement_edge_indices[level].end())
-    return true;
-  else
-    return false;
+  return !at_refinement_edge (level, index);
 }
 
 inline
@@ -274,9 +278,8 @@ MGConstrainedDoFs::at_refinement_edge (const unsigned int level,
                                        const types::global_dof_index index) const
 {
   AssertIndexRange(level, refinement_edge_indices.size());
-  AssertIndexRange(index, refinement_edge_indices[level].size());
 
-  return refinement_edge_indices[level][index];
+  return refinement_edge_indices[level].is_element(index);
 }
 
 
@@ -286,9 +289,8 @@ MGConstrainedDoFs::at_refinement_edge_boundary (const unsigned int level,
                                                 const types::global_dof_index index) const
 {
   AssertIndexRange(level, refinement_edge_boundary_indices.size());
-  AssertIndexRange(index, refinement_edge_boundary_indices[level].size());
 
-  return refinement_edge_boundary_indices[level][index];
+  return refinement_edge_boundary_indices[level].is_element(index);
 }
 
 inline
@@ -299,24 +301,54 @@ MGConstrainedDoFs::get_boundary_indices () const
 }
 
 inline
-const std::vector<std::set<types::global_dof_index> > &
-MGConstrainedDoFs::get_non_refinement_edge_indices () const
+const std::vector<std::vector<bool> > &
+MGConstrainedDoFs::get_refinement_edge_indices () const
 {
-  return non_refinement_edge_indices;
+  if (refinement_edge_indices_old.size()!=refinement_edge_indices.size())
+    {
+      unsigned int n_levels = refinement_edge_indices.size();
+      refinement_edge_indices_old.resize(n_levels);
+      for (unsigned int l=0;l<n_levels;++l)
+        {
+          refinement_edge_indices_old[l].resize(refinement_edge_indices[l].size(), false);
+          refinement_edge_indices[l].fill_binary_vector(refinement_edge_indices_old[l]);
+        }
+    }
+
+  return refinement_edge_indices_old;
 }
 
 inline
-const std::vector<std::vector<bool> > &
-MGConstrainedDoFs::get_refinement_edge_indices () const
+const IndexSet &
+MGConstrainedDoFs::get_refinement_edge_indices (unsigned int level) const
 {
-  return refinement_edge_indices;
+  AssertIndexRange(level, refinement_edge_indices.size());
+  return refinement_edge_indices[level];
 }
 
 inline
 const std::vector<std::vector<bool> > &
 MGConstrainedDoFs::get_refinement_edge_boundary_indices () const
 {
-  return refinement_edge_boundary_indices;
+  if (refinement_edge_boundary_indices_old.size()!=refinement_edge_boundary_indices.size())
+    {
+      unsigned int n_levels = refinement_edge_boundary_indices.size();
+      refinement_edge_boundary_indices_old.resize(n_levels);
+      for (unsigned int l=0;l<n_levels;++l)
+        {
+          refinement_edge_boundary_indices_old[l].resize(refinement_edge_boundary_indices[l].size(), false);
+          refinement_edge_boundary_indices[l].fill_binary_vector(refinement_edge_boundary_indices_old[l]);
+        }
+    }
+
+  return refinement_edge_boundary_indices_old;
+}
+
+inline
+const IndexSet &
+MGConstrainedDoFs::get_refinement_edge_boundary_indices (unsigned int level) const
+{
+  return refinement_edge_boundary_indices[level];
 }
 
 inline
@@ -328,20 +360,6 @@ MGConstrainedDoFs::set_boundary_values () const
   return boundary_values_need_to_be_set;
 }
 
-inline
-bool
-MGConstrainedDoFs::continuity_across_refinement_edges () const
-{
-  bool is_continuous = false;
-  for (unsigned int l=0; l<refinement_edge_indices.size(); ++l)
-    for (unsigned int i=0; i<refinement_edge_indices[l].size(); ++i)
-      if (refinement_edge_indices[l][i])
-        {
-          is_continuous = true;
-          return is_continuous;
-        }
-  return is_continuous;
-}
 
 DEAL_II_NAMESPACE_CLOSE
 
index e06cbc96d9cd8d07fa8ac6b3ff9da0336c49134c..bb2b9b08d156ce401bf98a315e52fc34945e97a2 100644 (file)
@@ -292,8 +292,8 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<std::vector<bool> >  &interface_dofs,
-                                std::vector<std::vector<bool> >  &boundary_interface_dofs);
+                                std::vector<IndexSet>  &interface_dofs,
+                                std::vector<IndexSet>  &boundary_interface_dofs);
 
   /**
    * Does the same as the function above,
@@ -302,7 +302,16 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<std::vector<bool> >  &interface_dofs);
+                                std::vector<IndexSet>  &interface_dofs);
+
+  /**
+   * As above but with a deprecated data structure. This makes one additional copy.
+   */
+  template <int dim, int spacedim>
+  void
+  extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
+                                std::vector<std::vector<bool> >  &interface_dofs) DEAL_II_DEPRECATED;
+
 
   template <int dim, int spacedim>
   void
index bc7ed01a807e4444b6f5398e927bef3fd6239800..bd85dd0cb0160a690ec54f293e1ac740abd35285 100644 (file)
@@ -410,7 +410,7 @@ namespace VectorTools
 
     const unsigned int max_rep_points = *std::max_element(n_rep_points.begin(),
                                                           n_rep_points.end());
-    std::vector<unsigned int>                 dofs_on_cell(fe.max_dofs_per_cell());
+    std::vector< types::global_dof_index>     dofs_on_cell(fe.max_dofs_per_cell());
     std::vector< Point<DH::space_dimension> > rep_points(max_rep_points);
 
     std::vector< std::vector<double> >           function_values_scalar(fe.size());
index 0536414ef5f7da6492603d0f392e533685b261c9..d07c32c912e3ca7fc17c47fdd9b6489e301d5467 100644 (file)
@@ -741,6 +741,27 @@ namespace TrilinosWrappers
 
 
 
+  BlockSparsityPattern::
+  BlockSparsityPattern (const std::vector<IndexSet> &row_parallel_partitioning,
+                        const std::vector<IndexSet> &col_parallel_partitioning,
+                        const std::vector<IndexSet> &writable_rows,
+                        const MPI_Comm              &communicator)
+    :
+    BlockSparsityPatternBase<SparsityPattern>
+    (row_parallel_partitioning.size(),
+     col_parallel_partitioning.size())
+  {
+    for (size_type i=0; i<row_parallel_partitioning.size(); ++i)
+      for (size_type j=0; j<col_parallel_partitioning.size(); ++j)
+        this->block(i,j).reinit(row_parallel_partitioning[i],
+                                col_parallel_partitioning[j],
+                                writable_rows[i],
+                                communicator);
+    this->collect_sizes();
+  }
+
+
+
   void
   BlockSparsityPattern::reinit (const std::vector<size_type> &row_block_sizes,
                                 const std::vector<size_type> &col_block_sizes)
index 4f21ce6136401bba2f25f5d8bba86e8c2fcde249..d1b5a6f5369e2b6428c6cec1f5413721f085cb49 100644 (file)
@@ -146,7 +146,8 @@ namespace TrilinosWrappers
     void
     BlockVector::reinit (const std::vector<IndexSet> &parallel_partitioning,
                          const std::vector<IndexSet> &ghost_values,
-                         const MPI_Comm              &communicator)
+                         const MPI_Comm              &communicator,
+                         const bool                   vector_writable)
     {
       const size_type no_blocks = parallel_partitioning.size();
       std::vector<size_type> block_sizes (no_blocks);
@@ -161,7 +162,8 @@ namespace TrilinosWrappers
         components.resize(n_blocks());
 
       for (size_type i=0; i<n_blocks(); ++i)
-        components[i].reinit(parallel_partitioning[i], ghost_values[i], communicator);
+        components[i].reinit(parallel_partitioning[i], ghost_values[i],
+                             communicator, vector_writable);
 
       collect_sizes();
     }
index 665c7ba1513884dfa8a454e23b2a50f38a0e406a..c486406f7fd7a92c1a65a1d3b08c8869521e3aaf 100644 (file)
@@ -389,6 +389,7 @@ namespace TrilinosWrappers
   void
   SparseMatrix::copy_from (const SparseMatrix &m)
   {
+    nonlocal_matrix.reset();
 
     // check whether we need to update the
     // partitioner or can just copy the data:
@@ -405,6 +406,9 @@ namespace TrilinosWrappers
         matrix.reset (new Epetra_FECrsMatrix(*m.matrix));
       }
 
+    if (m.nonlocal_matrix.get() != 0)
+      nonlocal_matrix.reset(new Epetra_CrsMatrix(Copy, m.nonlocal_matrix->Graph()));
+
     compress();
   }
 
@@ -475,6 +479,7 @@ namespace TrilinosWrappers
   {
     // release memory before reallocation
     matrix.reset();
+    nonlocal_matrix.reset();
 
     // if we want to exchange data, build a usual Trilinos sparsity pattern
     // and let that handle the exchange. otherwise, manually create a
@@ -582,9 +587,10 @@ namespace TrilinosWrappers
                   (Copy, sparsity_pattern.trilinos_sparsity_pattern(), false));
 
     if (sparsity_pattern.nonlocal_graph.get() != 0)
-      {
-        nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *sparsity_pattern.nonlocal_graph));
-      }
+      nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *sparsity_pattern.nonlocal_graph));
+    else
+      nonlocal_matrix.reset ();
+
     compress();
     last_action = Zero;
   }
@@ -598,6 +604,11 @@ namespace TrilinosWrappers
     matrix.reset ();
     matrix.reset (new Epetra_FECrsMatrix
                   (Copy, sparse_matrix.trilinos_sparsity_pattern(), false));
+    if (sparse_matrix.nonlocal_matrix != 0)
+      nonlocal_matrix.reset (new Epetra_CrsMatrix
+                             (Copy, sparse_matrix.nonlocal_matrix->Graph()));
+    else
+      nonlocal_matrix.reset();
 
     compress();
   }
@@ -755,6 +766,7 @@ namespace TrilinosWrappers
 
     const Epetra_CrsGraph *graph = &input_matrix.Graph();
 
+    nonlocal_matrix.reset();
     matrix.reset ();
     matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
 
@@ -801,8 +813,10 @@ namespace TrilinosWrappers
 
     // flush buffers
     int ierr;
-    if (nonlocal_matrix.get() != 0)
+    if (nonlocal_matrix.get() != 0 && mode == Add)
       {
+        // do only export in case of an add() operation, otherwise the owning
+        // processor must have set the correct entry
         nonlocal_matrix->FillComplete(*column_space_map, matrix->RowMap());
         Epetra_Export exporter(nonlocal_matrix->RowMap(), matrix->RowMap());
         ierr = matrix->Export(*nonlocal_matrix, exporter, mode);
index 85ca1111fcb3981a5710ef7cf35322dcec16e16e..b406c7300f20e1cc5b971e2e86a230337a80d345 100644 (file)
@@ -544,7 +544,7 @@ namespace TrilinosWrappers
     if (Utilities::MPI::n_mpi_processes(communicator) > 1)
       {
         Epetra_Map nonlocal_map =
-          nonlocal_partitioner.make_trilinos_map(communicator, false);
+          nonlocal_partitioner.make_trilinos_map(communicator, true);
         nonlocal_graph.reset(new Epetra_CrsGraph(Copy, nonlocal_map, 0));
       }
     else
index ea330e88e3a893acd1fb2254a17b809a41cc7ae5..4ff1d598e869a21466cbf2a8dbe67afba2af7ded 100644 (file)
@@ -186,6 +186,8 @@ namespace TrilinosWrappers
     Vector::reinit (const Epetra_Map &input_map,
                     const bool        fast)
     {
+      nonlocal_vector.reset();
+
       if (vector->Map().SameAs(input_map)==false)
         vector.reset (new Epetra_FEVector(input_map));
       else if (fast == false)
@@ -205,6 +207,8 @@ namespace TrilinosWrappers
                     const MPI_Comm &communicator,
                     const bool      fast)
     {
+      nonlocal_vector.reset();
+
       Epetra_Map map = parallel_partitioner.make_trilinos_map (communicator,
                                                                true);
       reinit (map, fast);
@@ -217,13 +221,11 @@ namespace TrilinosWrappers
                     const bool        fast,
                     const bool        allow_different_maps)
     {
-      // 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.
+      nonlocal_vector.reset();
+
+      // 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 (allow_different_maps == false)
         {
           if (vector->Map().SameAs(v.vector->Map()) == false)
@@ -250,14 +252,10 @@ namespace TrilinosWrappers
             }
         }
 
-      // Otherwise, we have to check
-      // that the two vectors are
-      // already of the same size,
-      // create an object for the data
-      // exchange and then insert all
-      // the data. The first assertion
-      // is only a check whether the
-      // user knows what she is doing.
+      // Otherwise, we have to check that the two vectors are already of the
+      // same size, create an object for the data exchange and then insert all
+      // the data. The first assertion is only a check whether the user knows
+      // what she is doing.
       else
         {
           Assert (fast == false,
@@ -284,19 +282,16 @@ namespace TrilinosWrappers
     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.
+      nonlocal_vector.reset();
+
+      // 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.
+      // create a vector that holds all the elements contained in the block
+      // vector. need to manually create an Epetra_Map.
       size_type n_elements = 0, added_elements = 0, block_offset = 0;
       for (size_type block=0; block<v.n_blocks(); ++block)
         n_elements += v.block(block).local_size();
@@ -349,26 +344,51 @@ namespace TrilinosWrappers
     }
 
 
-    void Vector::reinit(const IndexSet &local, const IndexSet &ghost, const MPI_Comm &communicator)
+    void Vector::reinit(const IndexSet &locally_owned_entries,
+                        const IndexSet &ghost_entries,
+                        const MPI_Comm &communicator,
+                        const bool      vector_writable)
     {
-      IndexSet parallel_partitioning = local;
-      parallel_partitioning.add_indices(ghost);
-      reinit(parallel_partitioning, communicator);
+      nonlocal_vector.reset();
+      if (vector_writable == false)
+        {
+          IndexSet parallel_partitioning = locally_owned_entries;
+          parallel_partitioning.add_indices(ghost_entries);
+          reinit(parallel_partitioning, communicator);
+        }
+      else
+        {
+          Epetra_Map map = locally_owned_entries.make_trilinos_map (communicator,
+                                                                    true);
+          Assert (map.IsOneToOne(),
+                  ExcMessage("A writable vector must not have ghost entries in "
+                             "its parallel partitioning"));
+          reinit (map);
+
+          IndexSet nonlocal_entries(ghost_entries);
+          nonlocal_entries.subtract_set(locally_owned_entries);
+          if (Utilities::MPI::n_mpi_processes(communicator) > 1)
+            {
+              Epetra_Map nonlocal_map =
+                nonlocal_entries.make_trilinos_map(communicator, true);
+              nonlocal_vector.reset(new Epetra_MultiVector(nonlocal_map, 1));
+            }
+        }
     }
 
 
     Vector &
     Vector::operator = (const Vector &v)
     {
-      // distinguish three cases. First case: both
-      // vectors have the same layout (just need to
-      // copy the local data, not reset the memory
-      // and the underlying Epetra_Map). The third
-      // case means that we have to rebuild the
-      // calling vector.
+      // distinguish three cases. First case: both vectors have the same
+      // layout (just need to copy the local data, not reset the memory and
+      // the underlying Epetra_Map). The third case means that we have to
+      // rebuild the calling vector.
       if (vector->Map().SameAs(v.vector->Map()))
         {
           *vector = *v.vector;
+          if (v.nonlocal_vector.get() != 0)
+            nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
           last_action = Zero;
         }
       // Second case: vectors have the same global
@@ -389,6 +409,9 @@ namespace TrilinosWrappers
           has_ghosts = v.has_ghosts;
         }
 
+      if (v.nonlocal_vector.get() != 0)
+        nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
+      
       return *this;
     }
 
@@ -397,6 +420,8 @@ namespace TrilinosWrappers
     Vector &
     Vector::operator = (const TrilinosWrappers::Vector &v)
     {
+      nonlocal_vector.reset();
+
       Assert (size() == v.size(), ExcDimensionMismatch(size(), v.size()));
 
       Epetra_Import data_exchange (vector->Map(), v.vector->Map());
index 48b0afcd0cef06b832a5e7062513aa82c5c6d615..95a9b1387baa949084d9cc0666fe08f2e8cf7c69 100644 (file)
@@ -21,6 +21,7 @@
 
 #  include <cmath>
 #  include <Epetra_Import.h>
+#  include <Epetra_Export.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -55,13 +56,9 @@ namespace TrilinosWrappers
       Assert (index < vector.size(),
               ExcIndexRange (index, 0, vector.size()));
 
-      // Trilinos allows for vectors
-      // to be referenced by the [] or
-      // () operators but only ()
-      // checks index bounds. We check
-      // these bounds by ourselves, so
-      // we can use []. Note that we
-      // can only get local values.
+      // Trilinos allows for vectors to be referenced by the [] or ()
+      // operators but only () checks index bounds. We check these bounds by
+      // ourselves, so we can use []. Note that we can only get local values.
 
       const TrilinosWrappers::types::int_type local_index =
         vector.vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
@@ -112,9 +109,8 @@ namespace TrilinosWrappers
   void
   VectorBase::clear ()
   {
-    // When we clear the vector,
-    // reset the pointer and generate
-    // an empty vector.
+    // When we clear the vector, reset the pointer and generate an empty
+    // vector.
 #ifdef DEAL_II_WITH_MPI
     Epetra_Map map (0, 0, Epetra_MpiComm(MPI_COMM_SELF));
 #else
@@ -186,18 +182,67 @@ namespace TrilinosWrappers
 
 
 
+  void
+  VectorBase::compress (::dealii::VectorOperation::values given_last_action)
+  {
+    //Select which mode to send to Trilinos. Note that we use last_action if
+    //available and ignore what the user tells us to detect wrongly mixed
+    //operations. Typically given_last_action is only used on machines that do
+    //not execute an operation (because they have no own cells for example).
+    Epetra_CombineMode mode = last_action;
+    if (last_action == Zero)
+      {
+        if (given_last_action==::dealii::VectorOperation::add)
+          mode = Add;
+        else if (given_last_action==::dealii::VectorOperation::insert)
+          mode = Insert;
+      }
+
+#ifdef DEBUG
+#  ifdef DEAL_II_WITH_MPI
+    // check that every process has decided to use the same mode. This will
+    // otherwise result in undefined behaviour in the call to
+    // GlobalAssemble().
+    double double_mode = mode;
+    Utilities::MPI::MinMaxAvg result
+      = Utilities::MPI::min_max_avg (double_mode,
+                                     dynamic_cast<const Epetra_MpiComm *>
+                                     (&vector_partitioner().Comm())->GetMpiComm());
+    Assert(result.max-result.min<1e-5,
+           ExcMessage ("Not all processors agree whether the last operation on "
+                       "this vector was an addition or a set operation. This will "
+                       "prevent the compress() operation from succeeding."));
+
+#  endif
+#endif
+
+    // Now pass over the information about what we did last to the vector.
+    int ierr = 0;
+    if (nonlocal_vector.get() == 0 || mode != Add)
+      ierr = vector->GlobalAssemble(mode);
+    else
+      {
+        Epetra_Export exporter(nonlocal_vector->Map(), vector->Map());
+        ierr = vector->Export(*nonlocal_vector, exporter, mode);
+        nonlocal_vector->PutScalar(0.);
+      }
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+    last_action = Zero;
+
+    compressed = true;
+  }
+
+
+
   TrilinosScalar
   VectorBase::el (const size_type index) const
   {
-    // Extract local indices in
-    // the vector.
+    // Extract local indices in the vector.
     TrilinosWrappers::types::int_type trilinos_i =
       vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
     TrilinosScalar value = 0.;
 
-    // If the element is not
-    // present on the current
-    // processor, we can't
+    // If the element is not present on the current processor, we can't
     // continue. Just print out 0.
 
     // TODO: Is this reasonable?
@@ -218,17 +263,13 @@ namespace TrilinosWrappers
   TrilinosScalar
   VectorBase::operator () (const size_type index) const
   {
-    // Extract local indices in
-    // the vector.
+    // Extract local indices in the vector.
     TrilinosWrappers::types::int_type trilinos_i =
       vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
     TrilinosScalar value = 0.;
 
-    // If the element is not present
-    // on the current processor, we
-    // can't continue. This is the
-    // main difference to the el()
-    // function.
+    // If the element is not present on the current processor, we can't
+    // continue. This is the main difference to the el() function.
     if (trilinos_i == -1 )
       {
         Assert (false, ExcAccessToNonlocalElement(index, local_range().first,
index a9da2451be3d470cb745a3b225f5ef32ed522042..85bf569a4fa13316f1800748a7ff0d3413e49f73 100644 (file)
@@ -1213,7 +1213,7 @@ namespace MGTools
         for (; cell!=endc; ++cell)
           {
             if (dof.get_tria().locally_owned_subdomain()!=numbers::invalid_subdomain_id
-                && cell->level_subdomain_id()!=dof.get_tria().locally_owned_subdomain())
+                && cell->level_subdomain_id()==numbers::artificial_subdomain_id)
               continue;
             const FiniteElement<dim> &fe = cell->get_fe();
             const unsigned int level = cell->level();
@@ -1247,7 +1247,7 @@ namespace MGTools
         endc = dof.end();
         for (; cell!=endc; ++cell)
           if (dof.get_tria().locally_owned_subdomain()==numbers::invalid_subdomain_id
-              || cell->level_subdomain_id()==dof.get_tria().locally_owned_subdomain())
+              || cell->level_subdomain_id()!=numbers::artificial_subdomain_id)
             for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
                  ++face_no)
               {
@@ -1416,7 +1416,7 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<std::vector<bool> >  &interface_dofs)
+                                std::vector<IndexSet>  &interface_dofs)
   {
     Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(),
             ExcDimensionMismatch (interface_dofs.size(),
@@ -1427,10 +1427,7 @@ namespace MGTools
         Assert (interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
                 ExcDimensionMismatch (interface_dofs[l].size(),
                                       mg_dof_handler.n_dofs(l)));
-
-        std::fill (interface_dofs[l].begin(),
-                   interface_dofs[l].end(),
-                   false);
+       interface_dofs[l].clear();
       }
 
     const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
@@ -1475,7 +1472,29 @@ namespace MGTools
 
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           if (cell_dofs[i])
-            interface_dofs[level][local_dof_indices[i]] = true;
+            interface_dofs[level].add_index(local_dof_indices[i]);
+      }
+  }
+
+  template <int dim, int spacedim>
+  void
+  extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
+                                std::vector<std::vector<bool> >  &interface_dofs)
+  {
+    std::vector<IndexSet> temp;
+    temp.resize(interface_dofs.size());
+    for (unsigned int l=0;l<interface_dofs.size();++l)
+      temp[l] = IndexSet(interface_dofs[l].size());
+
+    extract_inner_interface_dofs(mg_dof_handler, temp);
+
+    for (unsigned int l=0;l<interface_dofs.size();++l)
+      {
+        Assert (interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
+                ExcDimensionMismatch (interface_dofs[l].size(),
+                                      mg_dof_handler.n_dofs(l)));
+
+        temp[l].fill_binary_vector(interface_dofs[l]);
       }
   }
 
@@ -1483,7 +1502,7 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_non_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                              std::vector<std::set<types::global_dof_index> >  &non_interface_dofs)
+                              std::vector<std::set<unsigned int> >  &non_interface_dofs)
   {
     Assert (non_interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(),
             ExcDimensionMismatch (non_interface_dofs.size(),
@@ -1552,8 +1571,8 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<std::vector<bool> >  &interface_dofs,
-                                std::vector<std::vector<bool> >  &boundary_interface_dofs)
+                                std::vector<IndexSet>  &interface_dofs,
+                                std::vector<IndexSet>  &boundary_interface_dofs)
   {
     Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(),
             ExcDimensionMismatch (interface_dofs.size(),
@@ -1564,19 +1583,8 @@ namespace MGTools
 
     for (unsigned int l=0; l<mg_dof_handler.get_tria().n_global_levels(); ++l)
       {
-        Assert (interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
-                ExcDimensionMismatch (interface_dofs[l].size(),
-                                      mg_dof_handler.n_dofs(l)));
-        Assert (boundary_interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
-                ExcDimensionMismatch (boundary_interface_dofs[l].size(),
-                                      mg_dof_handler.n_dofs(l)));
-
-        std::fill (interface_dofs[l].begin(),
-                   interface_dofs[l].end(),
-                   false);
-        std::fill (boundary_interface_dofs[l].begin(),
-                   boundary_interface_dofs[l].end(),
-                   false);
+        interface_dofs[l].clear();
+        boundary_interface_dofs[l].clear();
       }
 
     const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
@@ -1639,10 +1647,10 @@ namespace MGTools
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           {
             if (cell_dofs[i])
-              interface_dofs[level][local_dof_indices[i]] = true;
+              interface_dofs[level].add_index(local_dof_indices[i]);
 
             if (boundary_cell_dofs[i])
-              boundary_interface_dofs[level][local_dof_indices[i]] = true;
+              boundary_interface_dofs[level].add_index(local_dof_indices[i]);
           }
       }
   }
index 8d4a68359b19f53818f037b3dadb6a6908689da3..7f25f462f0714b93e4086c04905e78466077b06e 100644 (file)
@@ -118,12 +118,17 @@ for (deal_II_dimension : DIMENSIONS)
     template
        void
        extract_inner_interface_dofs (const DoFHandler<deal_II_dimension> &mg_dof_handler,
-                                     std::vector<std::vector<bool> >  &interface_dofs,
-                                     std::vector<std::vector<bool> >  &boundary_interface_dofs);
+                                     std::vector<IndexSet>  &interface_dofs,
+                                     std::vector<IndexSet>  &boundary_interface_dofs);
+      template
+       void
+       extract_inner_interface_dofs (const DoFHandler<deal_II_dimension> &mg_dof_handler,
+                                     std::vector<IndexSet>  &interface_dofs);
       template
        void
        extract_inner_interface_dofs (const DoFHandler<deal_II_dimension> &mg_dof_handler,
                                      std::vector<std::vector<bool> >  &interface_dofs);
+
       template
         void
         extract_non_interface_dofs (const DoFHandler<deal_II_dimension> & mg_dof_handler,
index d8dcd64ed15f1e11ea67a7c999069a10492f9ad3..35442b79a0fcf656bf32de0f09954e7010b58089 100644 (file)
@@ -442,7 +442,7 @@ void MGTransferBlockBase::build_matrices (
                 constraints_per_block[block].resize(0);
                 constraints_per_block[block].resize(n_dofs, 0);
                 for (types::global_dof_index i=0; i<n_dofs; ++i, ++index)
-                  constraints_per_block[block][i] = constrain_indices[index] == 1;
+                  constraints_per_block[block][i] = (constrain_indices[index] == 1);
 
                 for (types::global_dof_index i=0; i<n_dofs; ++i)
                   {

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.