)
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)
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})
-
MARK_AS_ADVANCED(
atlas_LIBRARY
+ atlcblas_LIBRARY
+ atllapack_LIBRARY
blas_LIBRARY
eigen_blas_LIBRARY
+ f77blas_LIBRARY
gslcblas_LIBRARY
lapack_LIBRARY
m_LIBRARY
#
# 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)
# #
########################################################################
+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"):
#
# 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
* <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)
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
-// ---------------------------------------------------------------------
+ // ---------------------------------------------------------------------
// $Id$
//
// Copyright (C) 2013 by the deal.II authors
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>
}
+@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},
}
</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 -->
&AdvectionProblem<dim>::integrate_cell_term,
&AdvectionProblem<dim>::integrate_boundary_term,
&AdvectionProblem<dim>::integrate_face_term,
- assembler, true);
+ assembler);
}
BlockSparsityPattern (const std::vector<IndexSet> ¶llel_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.
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.
*/
}
-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
{
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;
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,
*
* @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;
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.
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.
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);
{
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;
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> ¶llel_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> ¶llel_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> ¶llel_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);
~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> ¶llel_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> ¶llel_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);
*
* @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;
/**
- * 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;
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);
}
+
inline
BlockVector::BlockVector (const std::vector<IndexSet> ¶llel_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)
{
}
+
inline
bool
BlockVector::has_ghost_elements() const
return ghosted;
}
+
+
inline
void
BlockVector::swap (BlockVector &v)
*
* 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
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]));
- // 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
* 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
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;
*/
//@{
/**
- * 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);
~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,
*/
//@{
/**
- * 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 ¶llel_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 ¶llel_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>
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 ¶llel_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 ¶llel_partitioning,
*/
//@{
/**
- * 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 ¶llel_partitioning,
const MPI_Comm &communicator = MPI_COMM_WORLD);
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 ¶llel_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 ¶llel_partitioning,
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 ¶llel_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);
//@}
};
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,
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;
};
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);
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;
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;
{
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;
//@{
/**
- * 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);
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.
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
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().
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;
//@}
//@{
/**
- * 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().
*/
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().
*/
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,
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().
*/
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);
//@{
/**
- * This collective set operation
- * is of lower level and can
- * handle anything else —
- * 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 — 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,
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,
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,
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,
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);
//@{
/**
- * 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,
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;
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;
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));
}
- 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 ()
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ if (nonlocal_vector.get() != 0)
+ nonlocal_vector->PutScalar(0.);
+
return *this;
}
{
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),
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)
}
}
{
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),
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;
+ }
}
}
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 ()
{
(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;
}
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;
}
const Worker &function;
const unsigned int partition;
const internal::MatrixFreeFunctions::TaskInfo &task_info;
+ const bool is_blocked;
};
} // end of namespace partition
};
- 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 ()
{
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;
}
const Worker &worker;
const unsigned int partition;
const internal::MatrixFreeFunctions::TaskInfo &task_info;
+ const bool is_blocked;
};
} // end of namespace color
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())
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;
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())
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
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;
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++)
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)
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
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++;
/// 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:
/**
* 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;
};
//----------------------------------------------------------------------//
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)
{
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)
{
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;
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)
{
* 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,
*/
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)
{
}
* 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;
- }
-
-
/**
* 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
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()
? 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
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);
// 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
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
}
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);
// 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);
}
* @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,
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);
}
/**
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>
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;
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);
+ }
}
{
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));
}
{
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
{
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);
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
{
// 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);
* 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:
/**
*/
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
*
* 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;
};
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);
}
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);
}
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();
}
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
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
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);
}
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
}
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
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
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,
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
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());
+ 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)
void
BlockVector::reinit (const std::vector<IndexSet> ¶llel_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);
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();
}
void
SparseMatrix::copy_from (const SparseMatrix &m)
{
+ nonlocal_matrix.reset();
// check whether we need to update the
// partitioner or can just copy the data:
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();
}
{
// 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
(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;
}
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();
}
const Epetra_CrsGraph *graph = &input_matrix.Graph();
+ nonlocal_matrix.reset();
matrix.reset ();
matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
// 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);
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
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)
const MPI_Comm &communicator,
const bool fast)
{
+ nonlocal_vector.reset();
+
Epetra_Map map = parallel_partitioner.make_trilinos_map (communicator,
true);
reinit (map, fast);
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)
}
}
- // 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,
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();
}
- 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
has_ghosts = v.has_ghosts;
}
+ if (v.nonlocal_vector.get() != 0)
+ nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
+
return *this;
}
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());
# include <cmath>
# include <Epetra_Import.h>
+# include <Epetra_Export.h>
DEAL_II_NAMESPACE_OPEN
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));
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
+ 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?
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,
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();
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)
{
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(),
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();
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]);
}
}
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(),
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(),
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();
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]);
}
}
}
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,
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)
{
projects / dealii-svn.git / commitdiff